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• ' Abstract: Matrix element reweighting is a powerful experimental technique widely em- 

f^ ' ployed to maximize the amount of information that can be extracted from a collider data 

^P ■ set. We present a procedure that allows to automatically evaluate the weights for any 

process of interest in the standard model and beyond. Given the initial, intermediate and 
final state particles, and the transfer functions for the final physics objects, such as leptons, 
jets, missing transverse energy, our algorithm creates a phase-space mapping designed to 
C^ \ efficiently perform the integration of the squared matrix element and the transfer functions. 

The implementation builds up on MadGraph, it is completely automatized and publicly 
available. A few sample applications are presented that show the capabilities of the code 
and illustrate the possibilities for new studies that such an approach opens up. 
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1. Introduction 

Along with the ongoing experimental activity in Run II at the Tevatron and at the now 
operational Large Hadron Collider at CERN, in the last years a significant effort by the 
high-energy community (including both theorists and experimentalists) has been devoted 
to devise new and more efficient strategies to identify physics beyond the standard model 
in collider data. In many of the new physics scenarios, new states exist at TeV scale that 
decay very quickly, and are not expected to leave any trace that can be reconstructed in 
the detector. Hence their existence and properties must be inferred from the distributions 
and properties of standard model particles that can originate from the decay of heavier 
not-yet-discovered resonances. 

The problem of identifying such decay patterns and from those of measuring the prop- 
erties of the new states is particularly intricate when the expected experimental signatures 
involve a complex final state, typically with several jets, leptons and missing energy. The 
latter, in particular, characterizes many of the models that aim at providing a candidate 
for dark matter consistent with the present observations. Several methods have been de- 
veloped during the last few years to improve on and eventually overcome this difficulty. 
For the purpose of identifying new physics, it is quite natural to first consider an approach 
that is not biased by strong theoretical assumptions, as the current knowledge of the vi- 
able theories that could lead to the production of new particles, is somehow limited. In 
this context, different methods have been proposed to measure the mass spectrum of the 
new states in a model-independent way: specific observables are suggested and built that 
are mostly sensitive to the masses of the new heavy resonances entering the decay chains. 
The final power of a given method is a balance between how well the information of the 
visible quantities is exploited to constrain the unknown masses and the dependence on the 
experimental and theoretical systematic uncertainties. Examples in the literature include 
the end-point method — that is based on the end-point regions of the invariant mass dis- 
tributions built from visible particles — , and the polynomial method — that attempts to 
reconstruct the whole event from the visible momenta — . These two methods can also be 
combined to give a better constraint on the mass spectrum. The Mt2 method [||, ^, and 
its generalizations/recent developments — though based on more complicated observables- 
also follow the same philosophy. For a recent review on these kinematics methods see ||3| 
and references therein. 

These model-independent methods will be determinant in constraining the mass spec- 
trum of new resonances. However, hy construction, most of them will not exploit or provide 
any information on other properties of the new particles, such as spin and coupling struc- 
ture. As another example, the precise measurement of the absolute mass of each particle 
entering a specific decay chain that ends with two missing particles remains challenging, 
especially in the case of short-length decay chains. In this context, it is useful to consider 
complementary and model-dependent tools for the investigation of properties of the new 
physics states. 

Matrix element reweighting [§, |5|, |6|, 0, P, ^, IC] is an example of such a method. It 



dramatically differs from the previous ones in that it requires at least one theoretical as- 



sumption as a starting point. Each assumption specifies the rules needed to compute the 
probability distribution associated with the process under study. Following a Bayesian ap- 
proach, the method assigns a probability to each hypothesis given a sample of experimental 
events, and in this way provides a discriminator among the different hypotheses. Another 
important feature of the method is that it makes maximal use of both experimental in- 
formation and the theoretical model (via the amplitude) on an event-by-event basis. This 
optimal use of the experimental data as well as the theoretical knowledge opens the door 
to new studies not only on mass measurements but also on the identification of the spin 
and coupling type of new particles. 

The matrix element method has been extensively exploited in the last decade. The 
best known example is its application to top-quark-pair production — investigated both 



by the CDF and D0 collaborations |11, 12, |l^, |lj, |T^ — which has led to the single most 



precise measurement of the top-quark mass. Recently, it has contributed to the observation 



of single-top production |16, 17] and to set an upper limit on the boson production of a 
standard model Higgs [|l8[, currently excluded in the mass region 158 GeV < mn < 175 
GeV at 95% confidence level ||]. 

In principle, the matrix element method is expected to provide a powerful investigation 
tool in many other analyses, in particular those dedicated to the search of new resonances 
and the study of their properties. However, in practice, its application is not straightfor- 
ward. In order to evaluate the weights to be attached to each experimental event, a difficult 
convolution of the theoretical information on the hard scattering (i.e., the matrix element 
squared) with the experimentally available information on the final state (encoded in the 
so-called transfer functions) has to be undertaken. The numerical efficiency (and therefore 
the speed) of such integration is currently a serious limitation. The matrix element squared 
as well as the transfer functions present variations by several orders of magnitude in differ- 
ent regions of the phase space. To overcome this difficulty, the integration technique has 
to be efficiently adapted to the shape of the integrand. To our knowledge, this problem 
has only been solved in very specific cases. 

In this work, we propose a general algorithm aimed at evaluating the weights appearing 
in the matrix element method. Given an arbitrary decay chain and the associated transfer 
function, our procedure first automatically assigns the optimized phase-space mappings 
designed to match as much as possible the peaks in the integrand, and then performs the 
phase-space integrations to evaluate the weights. Our implementation, which is fully auto- 
matic, is based on MadGraph, as it uses its matrix element amplitudes and the information 
on the topology of the diagrams. We dub the corresponding public code MadWeight. 

The paper is organized as follows. In Section |2| we review the basic features of the 
matrix element method. In Section ^ we expose our algorithm for the computation of the 
weights in the matrix element method. We present some illustrations in Section ^ and our 
conclusion in the last Section. 

2. The matrix element method 

As mentioned in the introduction, the matrix element method is a procedure to extract 



theoretical information, in the form of set of parameters a from a sample of experimen- 
tal events.^ Let us identify an event by the set x of experimentally available quantities 
(such as transverse momenta, rapidities, and so on). For each observed event a conditional 
probability P{x\Q.),i.e., a weight, is built that quantifies the "agreement" between the the- 
oretical framework a and the experimental event x. In the computation of the weights, one 
factorizes high-energy effects associated with the production of a parton-level configuration 
y into a calculable probability Pa{y). The evolution of the parton-level configuration y 
into a reconstructed event x in the detector is modeled by a transfer function W{x, y). As 
a result, the weight of a specific event x is of the form 

P{x\cx) = j dyP^{y)W{x,y). (2.1) 

In the specific case of a hadron collider, the parton-level probability Pa{y) can be expressed 
as a product of the squared matrix element |MQ,p(y), the parton distribution functions 
(pdf's) fi{qi) and f2{<l2) and the phase-space measure d^{y), such that the weight reads 

P(a,|a) = 1- [ d^{y)dqidq2h{qi)f2{q2)\M^\\y)W{x,y) . (2.2) 

0"a J 



The normalization by the total cross section a a in Eq. ( |2.2| ) ensures that P{x\oi.) is a 
probability density^: j P{x\cx)dx = 1. Once this probability density has been computed 
for each event Xi, the most probable value for a can be obtained through a likelihood 
maximization method. 

Eq. ( p.2| ) is central to this paper as it provides an explicit definition of the weight to be 
associated with a given event in terms of the convolution of tree- level matrix element, the 
pdf's and the transfer functions. One of the main working assumptions in the application 
of the matrix element method is that the transfer functions are " factorisable" ,i.e., they can 
be written as the product of single-particle resolution functions 

n 

W{x,y) = J{Wi{x\y'), (2.3) 

1=1 

where x* and y* stand for the measured quantities and the phase-space variables associated 
with the particle i, respectively. In practice, a further simplification is employed, where the 
transfer function associated with a single reconstructed object (such as a jet or a lepton) 
is written as a product of resolutions associated with the physical quantities measured in 
the detector: 

W,{x\y') = W^{x\y')Wl>{x\y')Wt{x\y'), (2.4) 



^Normally a. labels the different parameters in a given model (such as, for example, a mass or the value 
of a coupling). In this paper, however, we use a more general definition that also includes labelling different 
physics models. 

^We assume that the transfer function is also normalized to 1. 



where E, rj and (j) ^-re the reconstructed energy, rapidity and azimuthal angle, in most of 
the general purpose detectors the direction of a visible particle^ is well measured, so that 
the associated transfer function can be modeled by a narrow Gaussian. On the other hand, 
the resolution in energy strongly depends on the particle's type. For leptons it can be 
taken as a narrow Gaussian function whereas for jets a more involved parametrization of 
the resolution function is needed. 

In this work we provide a general solution to the problem of performing the integration 
in Eq. ( |2.2| ) in an efficient way. To better grasp the challenge that computing the integral 
in the numerator of Eq. ( |2.2| ) poses, it is useful to consider two limiting cases, where the 
problem simplifies. 

First let us imagine to have an "ideal" detector that could measure exactly the energies 
and momenta of all final state particles (including normally invisible ones), i.e., W{x, y) = 
6{x — t/).^ In this case no integration would be necessary and the weight in the numerator 
would be proportional to the corresponding squared matrix element, \Ma\'^{x). Nowadays, 
the determination of |Mo,p(a;) at the tree-level can be done automatically by several public 
codes and poses no difficulty. So apart from the normalization, discussed below, the weight 
calculation would therefore be trivial. 

As a second limiting case, one can also consider an ideal "no detector" option, i.e., 
choose the transfer function W{x,y) = 1. Then the integration would reduce to the 
computation of the total cross section, i.e., that of the denominator of Eq. (|2.2| ) as P{x\a) = 
1. This problem is not an easy one on its own: the matrix element has a very complicated 
peak structure, corresponding to the propagators of the Feynman diagrams being large. 
However, by observing that the leading peaks come from the sum of the squares of each 
diagram, together with the fact that it is always possible to find a parametrization of the 
phase space in terms of invariants that maps exactly those in the propagators |20|, makes 



the problem treatable (see for example Refs. |21, 22] and the discussion in the following 
Section). 

For a realistic detector the situation is in between the two above, where some par- 
ticles are well measured (charged leptons), other less (jets), and some completely missed 
(neutrinos). In this case the integration becomes extremely difficult as it involves an inte- 
grand with simultaneous peaks in sets of different variables that it is not possible, even in 
principle, to disentangle. 

3. Computation of the weights 

The evaluation of multi-dimensional integrals is often approached by standard adaptive 
Monte Carlo techniques. These techniques are well illustrated in the computation of total 
cross sections: phase-space mappings that "flatten" specific peaks in the integrand are 



^Throughout the paper, the expression visible particle refers to a lepton or a jet of which momentum is 
reconstructed in the detector. 

*For jets, the transfer functions also include genuine QCD effects like showering and hadronization. 
For the sake of the argument, we consider also "ideal" jets where the identification jet/parton is perfectly 
unambiguous. 



combined together in a multichannel integration. Here we also follow this approach. In our 
case, the phase-space mappings optimized for the computation of the weights in Eq. ( ^^ ) 
are rather involved because of the complex structure of peaks in the integrand, as it was 
discussed in the previous section. 

In this section we present our integration procedure and its implementation in a fully 
general algorithm. We first recall the basic principle of an adaptive Monte Carlo integration 



in Section 3.1 and then we describe the phase-space mappings optimized for the computa- 



tion of the weights in Section 3.2. We explain how we build a phase-space generator based 



on these new phase-space mappings in Section B^ and how we combine different phase- 



space mappings in a multi-channel integration in Section |3.4| . We validate our phase-space 
generator with several checks in Section |3.5|. 



3.1 Adaptive Monte Carlo techniques 

Adaptive Monte Carlo integration is a powerful numerical technique for the integration of 
a highly non-uniform function. It consists of sampling randomly the volume of integration 
according to a probability density that is adjusted iteratively to the shape of the integrand. 
The probability density is parametrized by a separable function 

p{z)=pi{z^)p2{z^) ...Pdiz'') (3.1) 

where each factor pi is a step function. If such a parametrization of the probability density 
function is appropriate to approximate the shape of the integrand, the adaptive integration 
procedure speeds up the convergence by increasing the density of evaluations in the regions 
where the integrand is large. In the case of a very sharp integrand, this condition is 
essentially fulfilled provided that the strength of each narrow peak in the integrand is 
associated with a single variable that in turn can be mapped onto one variable of integration 
z*. In that case, the integrand expressed in the parametrization z is of the form 




/(^) = M^') X R{^) (3.2) 



where the functions /j's may vary abruptly while the "remainder" non-factorisable function 
R{z) is essentially flat over the region under integration. 

If the integrand expressed in the phase-space mapping z presents a structure of sharp 
peaks that does not follow the factorized form in Eq. (p.2|), the adaptive integration pro- 
cedure is bound to fail. However, if enough information about the shape of the integrand 
is available, a first change of variables z — )• 2' = P{z) that rotates the axes of integration 
can sometimes be applied such that in the new phase-space mapping z' , the importance 
of each peak in the integrand is controlled by a single variable of integration. After this 
change of variables is applied, the integrand expressed in the new variables z' is of the form 
given by Eq. (|3.2|) , and the separable density function p{z') can be successfully adapted to 
the shape of the integrand. 

We will use the adaptive Monte Carlo integrator VEGAS [^] to carry out the in- 
tegration in Eq. ([2.21). Thus the efficiency in computing the weights will depend on the 



parametrization of the phase-space measure that is used in the adaptive Monte-Carlo in- 
tegration. The optimized phase-space mappings are such that for each narrow peak either 
in the transfer function or in the matrix element, the variable that controls the strength 
of that peak is mapped onto a single variable of integration in the parametrization of the 
phase-space measure, in which case the integrand expressed in that parametrization has 
the form given in Eq. (^Z 



3.2 The ne^v phase-space mappings 

For the computation of the weights, there is generally no simple phase-space parametriza- 
tion that maps all the peaks in the integrand and in which the boundaries of the phase- 
space volume can be easily expressed. Our strategy is to start from the following standard 
parametrization of the phase-space measure 

^,is [tJ \Pi\'^d\Pi\ sin 9ideid(l)i\ /o ^4x4 / , V^ I a o\ 

'^^=111 2E(2tt)^ j dqidq2i27r) 6 \pi+P2-2^Pj\, (3.3) 

where i = 3, . . .n labels the final particles. In this parametrization, the strength of each 
peak in the transfer function is already mapped onto a single variable of integration, whereas 
none of the propagator enhancement in the squared amplitude is. Identifying the Lorentz 
invariants associated with the Breit-Wigner resonances and expressing them as functions of 



the integration variables in Eq. (3^) is straightforward. The difficult task is then to invert 
these functions in order to derive a phase-space measure that is parametrized by both 
these Lorentz invariants and the variables mapping the peaks in the transfer function. 
Along with this inversion, the 6 function associated with energy-momentum conservation 



in Eq. (3^) has to be integrated out. The resulting phase-space mappings can then be 
used in an adaptive Monte Carlo integration to compute the weights. 

These optimized phase-space mappings can be defined by specifying the transformation 
of the phase-space measure parametrization in Eq. ( |3.3| ) from which they result. So in this 
Section, we will describe the expression of this transformation in a generic case, as it is a 
convenient way to introduce the new phase-space mappings. For an arbitrary process, the 
transformation that leads to the appropriate parametrization of the phase-space measure 
can be carried out by organizing the integration variables in the standard parametrization 



in Eq. (3^) into different subsets of variables to which a suitable change of variables is 
applied. Each subset of variables and its associated change of variables will be called a 
block in the following. 

The first phase-space block that needs to be identified is called the main block (MB), 
and it includes some of the integration variables appearing in Eq. ( |3.3| ) to which a trans- 
formation is applied so that the 6 function associated with energy-momentum conservation 
is integrated out. The same transformation may also map some invariants entering in the 
expression of specific propagators to new variables of integration in the expression of the 
phase-space measure. The identification of the main block and the form of the associated 
transformation of variables is discussed in the following Section. The integration variables 
appearing in Eq. (|3.3|) that do not belong to the main block also experience a transforma- 



tion that can be expressed in terms of secondary blocks, as explained in Section 3.2.2 
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Figure 1: Illustration of a decay chain with no missing particle: (a) the full topology, (b) the 
corresponding reduced diagram. The variables in the main block are written explicitly. 



3.2.1 Identification of the main block 

The main block includes a certain number of integration variables among those appearing 
in Eq. ( [3.31 ). These variables are adjusted as a function of all other kinematic quantities 
associated with the decay chain to enforce the conservation of energy and momentum. 

We start by discussing the choice of the main block in the case of two specific decay 
chains, and then generalize to the case of an arbitrary decay chain. We first consider 
a topology with no missing particle. An example of such a decay chain is illustrated in 
Figure |l|(a) . In that case, it is natural to include the initial proton momentum fractions of 
partons (called Bjorken fractions) in the main block, as the integrand does not show any 
sharp sensitivity in these variables. The angle of any visible particle should be excluded 
from the main block as it controls the strength of a narrow resolution function W^ or W^ 



in Eq. (2.4): it should therefore be maintained as an integration variable to ensure the 
integration of the associated peak. However, normally the resolution in energy is much 
poorer than the resolution in angles. For this reason, a relatively efficient choice is to 
complete the main block by adding two momentum variables \pj\ and |p,| of particles i 
and j that are relatively less constrained by the transfer function. In this example, the 
main block contains exactly four variables, and the effect of the variable transformation is 
to integrate out the 5 function with these four variables. 

We next move on to the case of a decay chain including missing particles in the final 
state. The phase-space variables associated with the momenta of missing particles are not 
directly constrained by the transfer function. Therefore they do not need to be mapped 
onto variables of integration in the phase-space mapping. We can identify the main block 
by selecting the momentum components of some missing particles instead of the energies 
of visible particles. A specific example of topology with missing particles is displayed in 
Figure §(a). One way to define the main block is to choose the set including the Bjorken 
fractions and the momentum components of the missing particle shown as a thick line in 



Px,Py,Pz 



(a) (b) 

Figure 2: Illustration of a decay chain with two missing particles (identified by the letter v): (a) 
the full topology, (b) the corresponding reduced diagram. The initial (resp. final) variables of the 
transformation associated with the main block are written explicitly, and the corresponding legs 
are shown as thick lines (resp. dashed lines). 

Figure 0(a). The change of variables associated with this main block remove these five 
variables from the set of integration variables in order to integrate out the 6 function in 
Eq. ( |3.3D and to map the invariant mass of the resonance decaying into the missing particle 
onto a variable of integration in the new parametrization of the phase-space measure. 

In order to generalize the discussion of the choice of the main block to the case of an 
arbitrary decay chain, it is useful to introduce the following representation of the main 
block and the corresponding transformation of variables: 

• In a branch of legs with no kinematic variable in the main block, the decay products 
of the initial particle in the branch are shrunk into a blob. 

• The variables in the main block are written explicitly and the corresponding legs are 
shown as thick lines. 

• The new integration variables resulting from the change of variables associated with 
the main block are also written explicitly, and the corresponding intermediate legs 
are shown as dashed lines. 

• All other intermediate legs that do not touch a blob are hidden behind a rectangular 
box. 

We refer to the resulting graph as the reduced diagram. As an illustration, the reduced 
diagrams for the two topologies shown in Figures |l](a) and Q(a) are displayed in Fig- 
ures |l](b) and |2|(b), respectively. In general a blob in a reduced diagram may hide a 
complicated branch of particles. But the change of variables associated with the main 
block is parametrized only by the total momentum of each blob, it does not depend on the 
structure inside the blob. 

We use these reduced diagrams to represent the main block in general. The minimum 
number of variables in the main block is four. After the variable transformation associated 
with the main block is applied, the 6 function in Eq. ( ^.3| ) is integrated out with these 
four variables that therefore do not appear in the new phase-space mapping resulting 



from this transformation. The main block may contain p > 4 integration variables. The 
transformation that is applied in that case removes all these p variables from the set of 
integration variables appearing in the parametrization of the phase-space measure, and 
introduces p — 4 new variables of integration. Each of these new variables map a Lorentz 
invariant that controls the strength of a specific propagator in the matrix element. Thus the 
variable transformation associated with the main block not only enforces the conservation 
of total energy and momentum, but also may possibly optimize the parametrization of the 
phase-space measure for the integration of some specific Breit-Wigner enhancements. 

As the variables in the main block are removed from the phase-space mapping af- 
ter the corresponding transformation is applied, an integration variable in the standard 
parametrization that controls the strength of a narrow peak in the transfer function is 
preferentially not included in the main block, otherwise the phase-space mapping after 
transformation would loose track of this variable and would be inappropriate for the inte- 
gration of the corresponding peak. From this observation, it is clear that the choice of the 
main block will act upon the efficiency of the Monte Carlo integration. 

Each of the main blocks treated in our code is illustrated by a reduced diagram in 
Figure ^. Their number is restricted because we only keep the main blocks for which the 
corresponding change of variables is invertible analytically. The corresponding formulas 
are discussed in Appendix ^. 

MB A. The transformation removes the Bjorken fractions qi and q2 and the norm of the 
three-momenta Pj , Pj of two visible particles from the set of integration variables in 
the parametrization of the phase-space measure. 
Example: pp — )• ZZ — )• 4j. 

MB B. The transformation removes the Bjorken fractions qi and g2 and the 3-momentum 
of a missing particle from the set of integration variables in the parametrization of 
the phase-space measure. The new integration variable is the invariant mass of the 
particle decaying into the missing particle. 
Example: pp — t- Z{W~^ — )• l^u). 

MB C. The transformation removes the Bjorken fractions qi and q2, the 3-momentum 
of a missing particle and the energy of a massless visible particle^ from the set of 
integration variables in the parametrization of the phase-space measure. The new 
integration variables are the Lorentz invariants m* and m* associated with the 
mother particles decaying into the missing and the massless particles, respectively. 
Example: pp — )• [i — )• h{W'^ — )• l~^v)]\i ^ b{W~ — )• jj)] with massless b quarks. 

MB D. The transformation removes the Bjorken fractions qi and 52 and the 3- momenta 
of two missing particles from the set of integration variables in the parametrization 
of the phase-space measure. The new integration variables are the Lorentz invariants 
m* , m* , m* and m* associated with the first and second mother particles of each 



''In the case the particle is massive, the corresponding change of variables turns out to be not analytically 
invertible. 
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Figure 3: The reduced diagrams representing the six main blocks that have been investigated in 
our procedure. 



missing particle. 
Example: pp — )• [i — )• b{W^ 



l+u) 



b{W- 



l-iy)]. 



MB E. The transformation removes the 3-momenta of two missing particles from the 
set of integration variables in the parametrization of the phase-space measure. The 
new integration variables are the Lorentz invariants m* and m* associated with the 
mother particles of each missing particle. The integration over the Bjorken fractions 
is expressed as an integration over the invariant mass and the rapidity of the colliding 
partons. 
Example: pp ^ H ^ {W+ -^ l+i^){W- -^I'v). 

MB F. The transformation removes the 3-momenta of two missing particles from the 
set of integration variables in the parametrization of the phase-space measure. The 
new integration variables are the Lorentz invariants m^ and rrv^ associated with the 
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Figure 4: A example of reduced diagram for wlricli the transformation that is applied to the MB 
cannot be inverted by means of analytical formulas. 



mother particles of each missing particle. 
Example: pp — )• (W^ — )• l~^i')(W~ — )• Z~z/). 

For a given decay chain, the main block can be chosen in several different ways, possibly 
with some more efficient than others. Roughly speaking a specific main block is appropri- 
ate for the computation of the weights provided that it does not contain a variable that 
controls the strength of a very sharp resolution function, and provided that no very sharp 
Breit-Wigner distribution is included in the square of the reduced diagram. This second 
condition comes from the fact that none of the invariants entering into the expression of 
the propagators inside the box is mapped onto a single variable of integration in the new 
parametrization of the phase-space measure, so that parametrization is not appropriate for 
the integration of the corresponding propagator enhancements. 

It should be stressed once again that each transformation of variables that is applied to 
the variables in main block has been implemented in the code analytically: for an arbitrary 
phase-space point, given the momenta of all the legs ending by the blobs and the invariant 
mass of each leg represented by a dashed line in the reduced diagram, the variables in the 
main block are determined by means of analytical expressions (see appendix |^) . We have 
not explored further the possibility to use numerical procedure for this step. So any change 
of variables that is not invertible analytically has been excluded in our algorithm. 

For example, the main block displayed in Figure ^ with three missing particles has not 
been considered. In principle, the 3-momenta of the three missing particles and the Bjorken 
fractions could be adjusted to satisfy eleven constraints induced by the seven resonances 
and the conservation of 4-momentum. As this adjustment cannot be done by means of 
analytical expressions, this case is dealt with the main block B or C in our procedure. 

3.2.2 Identification of the secondary blocks 

Once the main block has been defined and the corresponding transformation applied, the 
parametrization of the phase-space measure associated with the m external legs in all the 
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Figure 5: The four secondary blocks with their corresponding change of variables. Initial (resp. 
final) variables are written explicitly, and the corresponding legs are represented by thick lines (resp. 
dashed lines). 



blobs of the reduced diagram is still the standard one: 



n 



|Pjpd|pj| svnOidOi 
2E,{2TTf 



(3.4) 



The parametrization of the phase-space measure is further transformed by organizing the 
integration variables in Eq ( ^.4| ) into secondary blocks, i.e., into subsets of variables, each 
of them being subject to a specific change of variables. The change of variables of the 
simplest block is just the identity, in which case the variables in this block are maintained 
in the parametrization of the phase-space measure in Eq. (f3.4|). For the purpose of listing 
the other changes of variables that we have investigated, it is useful to represent a block 
and its corresponding change of variables by a diagram in the following way: 

• The variables involved in the transformation are written explicitly. The legs asso- 
ciated with the initial variables appear as thick lines. The legs associated with the 
final variables -which correspond to the invariants that enter into the expression of 
specific propagators- are shown as dashed lines. 

• A blob stands for a branch of legs of which total momentum parametrizes the change 
of variables related to the block. 

In this representation, a blob can a priori be itself decomposed into several secondary 
blocks. However, as in the case of the main block, the change of variables associated 
with a given secondary block is only parametrized by the total momentum of each branch 
represented by a blob, it does not depend on the details of these branches. The number of 
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Figure 6: Illustration of the structure in blocks optimizing the parametrization of the phase-space 
measure in the case of a specific decay chain. The missing particles are indicated by the Greek 
letter i/. 

implemented blocks in our algorithm is reduced by our requirement of analytically invertible 
changes of variables. These blocks are displayed in Figure |5[ The changes of variables 
associated with each block are discussed in Appendix pi 

SB A. The transformation removes the 3-momentum of a missing particle from the set 
of integration variables in the parametrization of the phase-space measure. The new 
integration variables are the Lorentz invariants m* 
first, second and third mother particles of this missing particle. 



m* and m* associated with the 



SB B. The transformation removes the energy and the polar angle of a missing particle 
from the set of integration variables in the parametrization of the phase-space mea- 
sure. The new integration variables are the Lorentz invariants m* and m* associated 
with the first and second mother particles of this missing particle. 

SB C/D. The transformation removes the energy of a missing particle from the set of 
integration variables in the parametrization of the phase-space measure (version C). 
The new integration variable is the Lorentz invariant m*^ associated with the mother 
particle of this missing particle. In version D of this block, the missing particle is 
replaced by a visible particle, but the transformation remains the same one. 

SB E. The transformation removes the momenta \pi\ and \p2\ of two visible particles pro- 
duced by the same resonance from the set of integration variables in the parametriza- 
tion of the phase-space measure. The new integration variables are the Lorentz in- 
variants 7T1,*, and m,*„ associated with the first and second mother particles of these 
visible particles. The corresponding change of variables is invertible analytically only 
if at least one of the two visible particles is massless. 

This completes the description of the blocks that can be used in our procedure to op- 
timize the parametrization of the phase-space measure for the computation of the weights. 
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As a example, we illustrate a composition in blocks in Figure y in the case of a specific 
decay chain. 

3.3 One-channel phase-space generator 

Given an optimized phase-space mapping defined by its structure in blocks, one can then 
consider a phase-space generator built upon this phase-space mapping. The generation of 
an arbitrary phase-space point proceeds in two steps: 1) the generation of the integration 
variables appearing in the optimized parametrization of the phase-space measure, 2) the 
determination of the momentum of each leg in the decay chain and the computation of the 
Jacobian factors. 

Concerning the first step, any variable of integration associated with the new phase- 
space mappings introduced in the previous section enters into one of the three following 
categories: 

1. The variable controls the strength of a resolution function. If the resolution function 
is a (5 distribution, the variable is fixed to the value associated with the experimental 
event. Otherwise, the grid of VEGAS is adapted such that the variable is generated 
according to a probability density that reproduces approximately the shape of the 
resolution function. 

2. The variable controls the strength of a propagator enhancement. In this case, the 
variable can be generated according to a probability density that reproduces exactly 
the shape of the propagator by using the inverse primitive function of a Breit-Wigner. 

3. The variable is either the polar or the azimuthal angle of a missing particle. In this 
case, the variable is generated according to a uniform distribution in the interval [0, vr] 
or [0, 27r] at the first iteration. The grid is adapted at each iteration to approximate 
the optimal probability density. 

Once the integration variables have been generated, the kinematics of the whole decay 
chain and the Jacobian factors are computed. For each block, the formulas that give 
the expression of the external momenta as a function of the variables of integration are 
discussed in the Appendix. These formulas are parametrized by the momentum of the 
branches represented by the blobs that appear in the graphical representation in Figures ^ 
and ^. For this reason, one needs to fill the kinematic variables in each block in a specific 
order, starting with the secondary blocks at the very end of the decay chain, and ending 
with the main block. 

This procedure is best illustrated with the example in Figure |6|. A phase-space point 
is defined by generating all the integration variables in the transformed expression of the 
phase-space measure: the invariant mass of each leg shown as a dashed line, the direction 
{0, (p) of any visible particle, and the energy of the visible particles represented by the solid 
thin lines. Then all other kinematic variables are determined as a function of the generated 
variables, first in the secondary blocks A and E, then in the secondary block D (by means 
of formulas that are parametrized by the kinematics in block E), and finally in the main 
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block B (by means of formulas that are parametrized by the kinematics of all the secondary 
blocks). Such an approach can be easily generalized to the case of an arbitrary decay chain. 
One difficulty in our approach is that the boundary of the physical phase-space vol- 
ume cannot be translated into simple conditions on the variables of integration. In order 
to implement the boundary, we simply check point-by-point in the phase-space that the 
variables in the main block are physical (for example, the Bjorken fractions qi or q2 can- 
not be larger than one), otherwise we throw away the phase-space point. In some cases, 
the fraction of unphysical points that are removed in this way may be large. Still, the 
algorithm is rather fast since the generation of a phase-space point is in general much less 
time-consuming than the evaluation of the squared matrix element. 

3.4 Mult i- channel phase-space generator 

If a given parametrization of the phase-space measure maps all the peaks in the integrand 
simultaneously, an adaptive Monte Carlo integration using only this channel is expected 
to be efficient. But most of the time, each peak in the integrand cannot be mapped onto a 
variable of integration in a single phase-space mapping, since the number of peaks is larger 
than the dimension of the phase space. ^ In these cases, we keep several channels, i.e. 
several phase-space parametrizations z ^ z' = Pi{z) such that each peak in the integrand 
is mapped onto a variable of integration in at least one channel. The total integration 
can be carried out using a multi-channel integration approach in which every channel i 
comes with a phase-space-dependent weight /3j(2;) > in the global parametrization of the 
phase-space measure: 

z^z' = P{z) = Y,fii{z)Pi{z), (3.5) 

i 

with the condition X^i ft(-^) ~ 1- Each weight I3i[z) must be chosen such that it is signifi- 
cant in the phase-space region where the corresponding channel Pi{z) is relevant. In the 
case of the computation of total cross sections, this condition can be automatically fulfilled 
by setting /3i{z) to be proportional to the amplitude squared of a single diagram associated 
with the channel i [£2|. In analogy to the single-diagram enhanced method, we choose to 
set the weight /3i{z) to be proportional to the product of the peaks that are mapped onto 
integration variables in the corresponding phase-space mapping. 

In comparison with previous implementations for the evaluation of the matrix element 
weights, this multi-channel approach is expected to speed up the convergence of the in- 
tegration, especially in the case of an over-constrained topology. An example of such a 



topology has been investigated in [12, 13, ^, where either the helicity of the W boson 
or the mass of the top quark is reconstructed from ti events in the semi-leptonic channel. 
In these analyses, a single channel was used for the evaluation of the weights, leaving un- 
mapped a subset of peaks in the integrand. On the contrary, our procedure always maps 
a given peak in the integrand onto a variable of integration in at least one channel. 



''This situation could also occur if one of the required blocks to build such a phase-space mapping 
corresponds to a change of variables that cannot be inverted analytically and hence has not been considered 
in our algorithm. 
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/ 


blocks 


integrated volume 


3 


MB A 


6.30 X 10-5 


3 


MB B 


6.30 X 10-5 


3 


MB C 


6.30 X 10-5 


6 


MB D 


694 GeV^ 


4 


MBE 


0.0166 GeV^ 


4 


MBF 


0.0166 GeV2 


5 


MB B + SB A 


3.89 GeV^ 


4 


MB B + SB B 


0.0166 GeV2 


3 


MB B + SB C 


6.30 X 10-5 


3 


MB B + SB D 


6.30 X 10-5 


4 


MB B + SB E 


0.0166 GeV^ 



Table 1: Phase-space volumes ^ dqidq2d(j)nl / {sqiq2) for I massless particles produced in hadron- 
hadron collisions at y^ = 1 TeV. The number I of final-state particles is indicated in the first 
column. The second and third columns indicate the structure in blocks defining the phase-space 
mapping that is used to calculate the volume with our phase-space generator, and the numerical 
value that we obtained. Each number is in agreement with the exact value of the phase-space 
volume at three digit accuracy. 

The whole procedure that we have presented so far has been implemented in the 
MadGraph framework, and the corresponding module has been named MadWeight. For a 
given decay chain and a transfer function for the final state objects, the optimized phase- 
space mappings are automatically selected, and the resulting multi-channel phase-space 
generator is used for the evaluation of the weights. While this procedure applies for virtually 
all cases, the speed of convergence of the numerical integration strongly depends on the 
process under investigation, and whether the calculation time is a serious limitation or not 
has to be assessed on a case- by-case basis. 

3.5 Validation of the phase-space generator 

One potential issue related to our phase-space mappings optimized for the computation of 
the weights is the fact that some of the associated Jacobians develop singularities in specific 
phase-space regions. These singular regions are an artefact of the change of variables. In 
our case they have a null measure in the integration volume. One can therefore split the 
integration volume into a volume Vi where the Jacobian is finite and a volume V2 that 
contains the singular region and that can be made arbitrary small compared to the volume 
Vi. At any given accuracy, we can ignore the contribution from the volume V2 provided 
that e = V2/V1 is sufficiently small. At the numerical level though, one may fear that 
instabilities will appear in this procedure. 

In practice, we have not encountered any numerical instabilities resulting from a change 
of variables that is associated with a specific phase-space block. Any phase-space block and 
the related change of variables that have been defined in our procedure have been checked 
by reproducing the volume of the entire phase-space region with our phase-space generator 
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I 


blocks 


integrated volume 


3 


MB A 


3.49 X 10-5 


3 


MB B 


3.49 X 10-5 


3* 


MB C 


4.13 X 10-5 


6 


MB D 


124 GeV^ 


4 


MBE 


8.17 X 10-3 GeV^ 


4 


MBF 


8.17 X 10-3 GeV^ 


5 


MB B + SB A 


1.28 GeV^ 


4 


MB B + SB B 


8.17 X 10-3 GeV^ 


3 


MB B + SB C 


3.49 X 10-5 


3 


MB B + SB D 


3.49 X 10-5 


4* 


MB B + SB E 


9.78 GeV^ 



Table 2: Phase-space volumes J dqidq2d(j)nl / {sqiq2) for I particles with a mass ?ti = 50 GeV 
produced in hadron-hadron collisions at y/s — 1 TeV. The number / of final-state particles is 
indicated in the first column. A star * indicates that the mass of one of the final state particles is 
set to zero, as this condition is required by one of the blocks. The second and third columns indicate 
the structure in blocks defining the phase-space mapping that is used to calculate the volume with 
our phase-space generator, and the numerical value that we obtained. Each number is in agreement 
with the exact value of the phase-space volume at three digit accuracy. 

using a parametrization of the phase-space measure that involves this block. This Monte 
Carlo procedure to compute the phase-space volume has a very poor convergence, as the 
phase-space mappings that are optimized for the computation of the weights are clearly 
inefficient for the computation of just the phase-space volume. Nevertheless, by increasing 
the number of generated phase-space points, we checked that the phase-space volume is 
reproduced with an accuracy better than one percent for each tested phase-space mapping. 
We first set the mass of the final-state particles to zero and obtained the results summarized 
in Table ||. We then considered the case of massive particles in the final state and obtained 
the results summarized in Table |2[ 

In order to validate the multichannel implementation, we also computed the total cross 
section of several processes by integrating the squared matrix element with our phase-space 
generator. This can be achieved by setting all transfer functions to one. Here again, the 
convergence of the numerical integration is poor, as the phase-space parametrization is not 
designed for such computation. By using a very high statistics, we reproduced the total 
cross sections associated with the processes listed in the first column of Table |3|. 



4. Example of applications 

In this section we illustrate a few examples of studies that can be achieved with Mad- 
Weight. The following analyses are based on simulated events generated with Mad- 
Graph/MadEvent (2! 



The events are passed through Pythia |26] for the showering and 
the hadronization. Electrons and muons are assumed to be reconstructed with 100% ef- 
ficiency and with excellent resolution if they have a pseudo-rapidity \r]\ < 2.4. Detector 
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process 


^MW/^ME 


channels 


blocks 


pp-^{W ^ jj)j 


0.982(6) 


3 


MBA 


pp-^ {W ^ Iv) 


0.9991(14) 


1 


MBB 


pp^ [W ^9r{f- >T-x)] 


1.003(5) 


1 


MB B; SB C 


pp -^ 2[fi -^ iix] 


1.020(5) 


3 


MB B,F; SB C 


pp -^ 2[t -^ b{W -^ li^i)] 


1.000(25) 


1 


MBD 


pp^[t^ b{W+ -^ lui)] + [i^ b{W- -^ jj)] 


0.94(5) 


6 


MB B, SB D,E 


pp^h^ {W+ -^ fl+Um){W' -^ 11- Urn) 


0.99(2) 


1 


MBE 



Table 3: Validation of the phase-space generator by computing total cross sections. The processes 
under consideration are written in the first column. The second column gives the ratio of the cross 
section computed with MadWeight over the one computed with MadEvent ||2^ . The third column 
indicates the number of channels that are used in the MadWeight integration, and the last column 
indicates the blocks that are involved in that integration. 



response simulation is performed using PGS |27] which takes into account geometrical 
acceptance, finite granularity and energy resolutions of typical calorimeters used in LHC 
experiments. Jets are then reconstructed based on the kt algorithm p 



P9| , pOl and applied 
on the calorimeter cells fired by the generated stable or quasi-stable particles. 

The transfer functions Wj (a;*,y*), Eq. ( p.4D with i running over all reconstructed jets 
are determined from an independent tt sample where well separated jets (including light and 
b jets) are matched to the corresponding partons. We consider a double- Gaussian shape 
function characterized by 5 parameters: the means and the widths of the two Gaussian 
distributions, and their relative normalization. We fit these five parameters in each 20 
GeV bin in jet energy from 40 GeV to 200 GeV. The energy dependence of the mean 
and the width of each Gaussian distribution is then approximated by the parametrization 
Cl+C2^fE+c^E, with the coefficients ci, C2 and C3 extracted from a x^ ^t to the values of the 
four parameters of Gaussian distributions in each energy bin. The relative normalization 
of the two Gaussian distributions is assumed to be energy independent, and is fixed to 
the average of the corresponding values in each energy bin. The typical resolution for jet 
energy is between 5 and 12 GeV, with tails parametrized by Gaussian of variances as large 
as 30 GeV. 



4.1 Top-quark mass measurement 

The top-quark mass measurement by means of the matrix element method was published 
for the first time by the D0 collaboration using the single-leptonic final state arising from 



top-quark pair production |11]. The method has been later extended also to include a 
simultaneous determination of the Jet-Energy-Scale uncertainty. The accuracy of the ex- 
perimental determination of mt has been further improved by the contribution of other 
studies based on matrix element method [^, 14, |l^. In these analyses, a dedicated phase- 
space integration was performed to define the event weight. Our algorithm provides this 
weight automatically based on the blocks given in Table |3|, last but one line. 

As an example of application of our automatic reweighting algorithm, we illustrate the 
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Figure 7: (a) Logarithmic likelihood values for a sample of 20 events generated with rrit^ input = 170 
GeV. The solid line is a parabolic fit to the points near the minimum. The statistic error is 
estimated by the half width of the distribution at log(L/imax) = 0.5 and is extracted from the 
fit. (b) Calibration of the matrix element mass fitting procedure. The errorbars correspond to the 
the value of mass of the top quark and the associated statistic error reconstructed from tt samples 
generated with different input values of rrit- The solid line is a linear fit to the four points and the 
dotted line corresponds to mice = f^input- 

performance of the method for the determination of the top-quark mass at the LHC, by 
using a small statistics of tt events in the single lepton final state: 



p'p 



[i^ b{W- ^ t,-u^][t ^ biW+ ^ jj)]. 



(4.1) 



For the sake of simplicity, we assume that there is no background and 20 signal events 
after selection. Pseudo-data have been simulated with an input top-quark mass at 170 
GeV. The selection requires one muon with a reconstructed transverse momentum above 
10 GeV and exactly four isolated jets with a reconstructed transverse momentum above 20 
GeV. 

The determination of the top-quark mass from our sample of pseudo-data is obtained 
by the minimization of — log(L) with respect of rrit where the likelihood L is defined -up 
to a normalization factor- by the product of the weights calculated for each event 



N 



log(^) = -^ log[P{xi;mt)] . 



(4.2) 



i=l 



The acceptance of the detector and the cuts imposed on the sample may depend on the 
input mass of the top quark. Such a dependence might introduce a bias in the extraction 
of m-top from the fit of the likelihood given in Eq. ([4.2|). We explicitly tested that in our 
pseudo-data this bias is very small and we therefore ignored it in this example. 

The values of — log[L(?ni)] for different assumptions of mt are displayed in Figure |^( a). 
A clear minimum is observed close to the input mass value. A parabolic fit gives the 
value rrit = 171.9 ± 2.0stat GeV. Using ten independent samples generated under the same 
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conditions, we extracted the expected mt value and the expected statistic error for such a 
measurement: mt = 173.5 ±3.7 stat GeV. We also identified a small significant shift between 
the input and the reconstructed mass by repeating the analysis for other input masses, as 
can be seen from the calibration curve displayed in Figure |^(b) . This curve was obtained by 
generating tt samples with top quark masses of 160, 170, 180 and 190 GeV with the same 
selection procedure and the same fitting procedure, and with a statistics of 100 events per 
sample. The small bias between the input and the reconstructed mass of the top quark can 
possibly result from the effect of initial state radiation, which is not completely removed 



in our event selection criteria |31]. 



4.2 Spin identification in decay chains with missing energy 

As a second illustration we address the challenge of determining the spin of new particles. 
A simple example is the production of a light charged Higgs boson with a mass close to the 
W boson mass. The information on the spin of a resonance is passed through the angular 
distribution of its decay products. If the momentum of each final-state particle produced 
in the decay chain is measured, the angular distributions can be reconstructed and the 
spin of the resonance identified. In fact in many cases, such as the one we have chosen, 
the final state is characterized by missing transverse energy from undetected particles 
and the angular distributions of the decay products cannot be fully determined. The 
interesting question becomes therefore whether the available information from the final 
state is sufficient to discriminate between different spin assignments in the decay chain. 
The matrix element method appears to be particularly relevant in this case, since the event 
weight will encompass the whole available event kinematics including the spin correlation 
effects that survive after the experimental reconstruction of the events. 

In the following example, we assume that the production of the signal and its irreducible 
background proceeds exclusively via the production of top quark pair that subsequently 
decay into H~^ + 6 or into W + b, with mfj± ~ m^r± , followed by a leptonic decay of both 
bosons. The signal process is 

pp^[t^ b{H+ -^ T+Vr)\[i^ b{W- -^ ^l-u^,)l (4.3) 

and the corresponding irreducible background results from the production of a pair of W 
bosons 

pp^[t^ b{W^ -^ T+Ur)][i^ b{W- -^ ^i'i>^,)]. (4.4) 

At the reconstruction level, we required the presence of exactly two jets with a pT larger 
than 20 GeV, one r"*" — assumed to be reconstructed as precisely as the other charged 
leptons — and one /x~. The cuts on these leptons are |ry| < 2.4 and px > 5 GeV. We 
reject the events containing photons or electrons with pj- > 5 GeV and |?7| < 2.4, while 
no restriction is imposed on the number of jets with a px less then 20 GeV. The transfer 
functions associated with the other particles have the same parametrization as described in 



Section 4.1. We arbitrarily choose a final relative normalization of signal and background 



events, working with a sample of 240 signal events and 760 background events. 
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Figure 8: Expected normalized distribution of events with respect to the discriminant d built upon 
(a) the matrix element weight and (b) the px of the tau for a pure signal sample (solid histogram) 
and for a pure background sample (dashed histogram). The errorbars are the distribution associated 
with the pseudo-experiment sample, assuming that the statistical error on the number N of events 
in a bin is given by VN. 



Before proceeding further, we stress that while providing an interesting case study, our 
example cannot be regarded very realistic. First, the relative cross sections and reconstruc- 
tion efficiencies for the signal and the background have been chosen arbitrarily. Second, 
this such a light charged Higgs is not favoured by the present constraints, which point to 
much higher masses. Finally, the tau lepton reconstruction is idealized as it is considered 
here on the same footing as as a muon. A more realistic approach would consist of taking 
into account the energy loss from the tau decay with a dedicated transfer function for the 
energy of the tau. Nonetheless, as shown below, this example illustrates quite well the 
power of the matrix element method. 

Let us define Psix), Pb{x) as the weights evaluated for the event final state x under 
the signal and the background hypotheses, respectively. These weights can be calculated 
from the signal and background full matrix element as defined in Eq. ( |2.2| ). Alternatively, 
they can be associated with a normalized differential cross section with respect to a single 
observable, such as the r"*" transverse momentum 



Ps,b{x) 



1 das,B 



[PTir"-)] 



(4.5) 



(^S,B dpT 

which also captures the spin effects. The advantage of the weights defined in Eq. (|4.5|) is 



in their simplicity: their evaluation only requires to use a standard phase-space generator 
that is optimized for the computation of cross sections. Such an observable, for example, 
is very commonly used in the determination of the polarization of the W bosons in top 
events and provides us with a useful benchmark to study the increased sensitivity that the 
matrix element method might provide. 

These weights can then be combined to build an event-by-event discriminating variable 

Ps{x) 



d{x) 



Ps{x) + Pb{x) 



(4.6) 
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Figure 9: x^ values associated with the fit of the pseudo-experiment data to the theoretical 
prediction parametrized by the fraction r of signal events for the weights calculated from (a) the 
matrix element, (b) the pt of the tau. Each dashed line represents the one-standard-deviation 
interval defined by the condition x^(r) < x^j^^ -I- 1. 



The normalized-to-one distributions of events as a function of the discriminant variable d 
are shown in Figure ^ for the two cases. The solid (resp. dashed) histogram is the distri- 
bution expected for a pure sample of signal (resp. background) events. These distributions 
have been generated from large samples, in order to allow us to neglect the statistical 
fluctuations. Spin correlation effects are expected to give rise to different values for the 
weights under the two spin hypotheses. Nevertheless, for most of the events, this disparity 
is expected to be small, resulting in a discriminant close to d ~ 0.5. Note that in our 
example, signal and background events are characterized by the same topology with inter- 
mediate particles of the same mass. Only the spin of the intermediate W or H resonances 
differ between the two decay chains. The distributions clearly shows that the discriminant 
power is substantially reduced when only the information on the transverse momentum of 
the T~^ is retained. 

Yet, for some events the discriminant is significantly different of 0.5, corresponding to 
configurations clearly favoured by one of the two hypotheses. Such events influence the 
shape of the distributions and allow us to distinguish them. One can take advantage of this 
difference to find out the fraction of signal events in the pseudo-experiment sample. The 
normalized-to-one distribution associated with the pseudo-experiment sample as a function 
of the discriminant variable d is also displayed in Figure ^. The fraction of signal events in 
the pseudo-experiment sample can be reconstructed by a least-square fit, i.e. by minimizing 



x'(r) 



bins in d 



{-Pdatajd) - [rVsjd) + (1 - r)VBid)]y 

[AVdataidW 



(4.7) 



where Vs, Vb, are the expected binned distributions for signal and background events and 
Vdata is the binned distribution associated with the pseudo-experiment sample. 

The x^ values as a function of the fraction of signal events are shown in Figure p. The 
best fit for each discriminant are obtained for r = 24 it 9% and r = 30 it 23%, respectively. 
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Figure 10: Generic decay chain corresponding to the production of two resonances followed by 
their decay into weakly-interacting and standard model particles. The weakly-interacting particles 
are represented by the dashed lines. 

Both results are compatible with the true fraction of signal events, but the discriminant 
using matrix elements increases the accuracy by more than a factor of two. 

4.3 Smuon pair production at the LHC 

Over the past fifteen years, a tremendous amount of work has been devoted to new tech- 
niques for mass reconstruction of new particles that might be produced at the LHC. Accord- 
ing to most scenarios, the hypothetical new physics states are not expected to be directly 
observed experimentally, i.e., they appear as intermediate states in specific decay chains or 
they escape from the detector without interacting with it. Their mass can hence only be 
reconstructed indirectly, by making a number of assumptions on the decay chain at work. 
The number of assumptions in turn is directly correlated to the amount of information 
that can be extracted from the decay chain. However, due to the lack of constraints on 
physics beyond the standard model, the proposed techniques have to be general enough, 
at least if they are aimed at reconstructing the mass of new hypothetical particles in the 
early stages of investigation. Further more, the limited knowledge of the detector has to 
be taken into account. It is in this context that a number of mass measurement techniques 
based on kinematic methods have been proposed in the literature. They can be classified 
according to the type of decay chains that they address and according to the assumptions 
on which they rely p]. 

Despite the plurality of kinematic variables that have been proposed, mass determi- 
nation remains very challenging for specific decay chains. One well-known example of a 
difficult topology is the production of two resonances followed by their decay into a weakly- 
interacting and a standard model particles, shown in Figure |l^. In this case, the kinematic 
methods that have been proposed to reconstruct simultaneously the mass of the two new 
particles require a very high statistics. Whether their sensitivity is sufficient under real 
experimental conditions still remains to be determined. A complementary way to address 
the same problem is to ask what would be the maximum sensitivity achievable, given a 
very detailed set of hypothesis to be tested that not only include masses but also the spin 
and couplings information, i.e., taking into account the full theoretical model prediction. 

The problem can be investigated with the matrix element method that usually makes 



use of the strongest assumptions on the analysed events [32|. One way to dramatically 



increase the theoretical information is to assume that the masses of the new physics states 
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Figure 11: Logarithmic likelihood as a function of the hypothesis values for [(tti- — 
TO? )/2'mjx^,my--^] built upon (a) the matrix element weights, (b) the transverse momentum of 
the fi^ and the invariant mass of the muons. 

are the only unknown properties of the decay chain. We therefore consider a specific decay 
chain corresponding to the topology in Figure |l^: the production of a pair of smuons 
followed by their decay into a muon and a neutralino 

pp^ (/i+ ^^+Xi)(Ar ^M^Xi)- (4.8) 

We suppose that we have isolated a pure sample of events that correspond to the decay 



chain in Eq. (4.8). Further, we assume a perfect reconstruction of the kinematics of the 
two muons in each event. Within these assumptions, the significance that can be achieved 
with the matrix element method provides us with an upper bound on the significance that 
can be delivered by any realistic analyses at a given luminosity. 

We have considered the following input hypothesis for the masses of the sparticles: 

m-/i,,mput = 150 GeV, m^,, input = 100 GeV. (4.9) 

Under this hypothesis, we have generated events corresponding to the decay chain in 
Eq. ([4.^). We have built a sample of fifty events with exactly one fi'^ and one ^~ with a 
transverse momentum larger than 20 GeV, and no other particles except maybe some jets 
with a pt less than 20 GeV. These events are regarded as a pseudo-experiment sample in 
the following. 

The sensitivity that can be achieved with the matrix element method has been anal- 
ysed by computing the weights P(xj|m^^,?TT,^^) for each event Xi in the pseudo-experiment 
sample. A bias may be introduced by the acceptance cuts. This effect has been corrected 
by normalizing the probability density in the acceptance region. This amounts to replace 
the factor 1/aa by the factor l/a^^ in the definition of the probability density in Eq (^), 
with a the cross section in the acceptance region. In terms of the weights normalized in 
the acceptance region, the unbiased likelihood has the usual form 

iV=50 



log L{mf,^ 



y^ logP(xj|mp 



(4.10) 
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It is advantageous to express the likelihood in term of the variable (m? — m3- )/2m^^, 
that corresponds to the momentum of each final state particle in the rest frame of the 
smuon from which it originates. The complementary variable can be chosen to be "i^^. 
The likelihood for different theoretical hypotheses is shown in Figure |ll|(a). The optimal 
value for the variable (fn? — m|^)/2m,^^ is 42 GeV, which corresponds to the input value. 
There is a very mild sensitivity with respect to variation of the complementary variable. 

One way to highlight the increase of sensitivity by using the complete theoretical and 
experimental information is to compare the profile of the likelihood built upon the matrix 
element weights (shown in Figure 0(a)) with the likelihood profile that is obtained by 
keeping only the information contained in the transverse momentum px^ of the /i"*" and 
the invariant mass M^^ of the muons. In order to simplify the computation of this second 
likelihood profile we neglect the correlations between the two variables pxfi and M^^. Thus 
the weight attached to each event is reduced to 

, I , 1 da . . \ 1 '^'^ / 71 J- I \ /A -,-,\ 

The resulting likelihood profile is displayed in Figure pT|(b). The comparison with Fig- 
ure |ll|(a) shows that the sensitivity to the theoretical hypothesis (m^^ , rrij^^ ) is dramatically 
reduced when only the information contained in the kinematic variables pT^ and M^^ is 
used. 

5. Conclusion 

We have presented a new algorithm that allows the automatic computation of the weights 
appearing in the matrix element method. Given an arbitrary decay chain and a transfer 
function tuned to the resolution of the detector, our code produces a specific phase-space 
generator that combines different phase-space mappings optimized for the integration of 
the product of the matrix element and the transfer function. The mappings are obtained 
by applying a specific transformation to the standard parametrization of the phase-space 
measure. This transformation is expressed in terms of a composition of changes of variables 
acting on different kinematic sectors of the topology. As a result, our algorithm leads to 
a modular structure, and hence is particularly convenient for future improvements. For 
example, the current implementation could be easily extended to include non-analytical 
changes of variables. 

The availability of a tool such as MadWeight that provides the resource for the auto- 
matic evaluation of the weights, sparing the user to focus on the technical details of matrix 
element generation and integration over phase space, paves the way to a potentially large 
number of new applications. First MadWeight could be used to improve our understanding 
of the matrix element method itself and of its limitations. For example, its implementation 
within Madgraph provides all the required computational tools to analyse the infiuence 
of additional jet radiation in a specific measurement, or to estimate the systematic error 
resulting from the parametrization of the transfer function. Second not only measurement 
of masses or cross sections in Standard Model could be achieved in a effortless and more 
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efficient way, but the matrix element method could also be employed in the search and 
identification of new physics models. Different hypotheses, such as those corresponding to 
(any) new physics scenario and/or benchmark parameter points could be tested against 
data and be assigned a meaningful relative probability. We look forward to exciting new 
developments in these directions. 
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A. Phase-space measure associated with the main blocks 

A.l MB A 

The notation for the phase-space variables associated with this main block is given in 



Figure 12. The two momenta pi and p2 correspond to the visible particles that enter into 

'?! -^ / / bil 
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|P2l 



Figure 12: Notation for the kinematics of MB A. 

the main block, along with the Bjorken fractions qi and 52- The standard phase-space 
parametrization associated with this main block reads 

The four-vectors Pjn and Pgn refer to the total momenta in the initial and final states, 
respectively. In our procedure we apply a change of variables that leads to the following 
parametrization of the phase-space measure 

167r2^ E dMhde2d^2 x J, (A.2) 
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where 9i and (j)i refer to the polar and azimuthal angles of particle i with respect to the 
beam axis. The Jacobian J of this transformation reads 

2, 



J 



iPil IP2I I cos (/)i sin 02 



sm (pi cos (p2 



(A.3) 



where s is the squared invariant mass of the colliding hadrons. The energies Ei, E2 of 
the final particles in the main block are adjusted to balance the transverse momentum 



^branches q£ ^^^j ^j^g branches represented by the blobs in Figure 12. We assume that this 



transverse momentum is different from zero, except maybe in a region of null measure (see 
the discussion in Section |3^ ). This requires the number of particles in the final states to be 
larger than or equal to three. Then the variables \pi\, IP2I can be expressed as the solution 
of the following linear system 



|p;^| sin^i cos 01 + IP2I sin^: 
\pi\ sin 01 sin 01 + ' 



'2 cos I 



^branches 
'2 = -Px 



|p2|sin02sin02 = -p!;'"'"'''''^ 



(A.4a) 
(A.4b) 



The Bjorken fractions (71,(72 are then fixed by imposing the conservation of total energy 
and total momentum along the beam axis. 

A.2 MB B 

The notation for the phase-space variables associated with this main block is given in 



Figure 13. The momentum pi corresponds to the missing particle that belongs to the main 
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Figure 13: Notation for the kinematics of MB B. 

block, along with the Bjorken fractions qi and g2- The momentum p2 corresponds to the 
branch that is directly connected to that missing particle. The variable S12 is the invariant 
(pi + ^2)^- The standard phase-space parametrization associated with this constrained 
sector reads 



dqidq2 



-{2nr6^H^-Pi 



fin J 



(A.5) 



(27r)32^i 

The four-vectors Pjn and Pgn refer to the total momenta in the initial and final states, 
respectively. In our procedure we apply a change of variables that leads to the following 
parametrization of the phase-space measure 



4ttEi 



dsi2 X J. 



(A.6) 
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The Jacobian J of this transformation is given by 



J = 1^22-^1 - E2PIZ 



(A.7) 



where s is the squared invariant mass of the colhding hadrons. The transverse momentum 
the missing particle is fixed by requiring that it balances the transverse momentum of all 



the branches represented by the blobs in Figure |13|. The component piz of momentum 
along the beam axis is fixed by imposing the invariant mass condition 

{Pl+P2f 



Sl2- 



(A. 



If the energy Ei of the missing particle is treated as an independent parameter, the left side 
of Eq. (|A.8| ) is a first-order polynomial in piz . We therefore obtain a unique expression for 
piz in terms of Ei. The mass-shell condition associated with the missing particle gives rise 
to up to two solutions for the energy Ei. Each solution that gives a real positive value for 
El and that leads to values of the Bjorken fractions qi and §2 between and 1 is kept, as 



it corresponds to a distinct physical phase-space point at which the Jacobian in Eq. ([A.?]) 
and the integrand must be evaluated. 

A.3 MB C 

The notation for the phase-space variables associated with this main blob is given in Fig- 



ure 



14. The momentum of the missing particle is denoted by pi, the momentum of the 
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Figure 14: Notation for the kinematics of MB C. 

branch directly connected to the missing particle is denoted by p2^ the momentum of the 
massless visible particle in the main block is denoted by ^3. The Bjorken fractions are 
denoted by qi and 52- The variables S12 and S123 refer to the invariants {pi + ^2)^ and 
[pi + P2 +^3)^1 respectively. The standard phase-space parametrization associated with 
this main block reads 

d^pi d?Pz 



dqidq2 



(27r)^5* (Pin - Pfin) 



(A.9) 



' (27r)32£;i i2TT)^2E3 

The four-vectors Pjn and Pgn refer to the total momenta in the initial and final states, 
respectively. In our procedure we apply a change of variables that leads to the following 
parametrization of the phase-space measure 



16tt'^EiE3 



id9sdsi2dsi23 X J- 



(A.IO) 
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The Jacobian J of this transformation is given by 



J = sin03 



^l^i 



XE2PIZ - XElP2z 



s 

2 cos((/)3) 003(63) E2P1XE3 sin(6'3) + 2 cos(03) cos{93)Eip2xEz sin(6'3) - 
2cos{<i)3)pizP2xE-i sin(6'3) + 2 (zos{4>3)pixP2zE3 sin(6'3) - 
2 cos{e3)E2PiyE3 sm{(j)3) sin(6'3) + 2 €08(63) Eip2yE3 sm{(f)3) sm{63) - 
2pizP2yE3 sm{(j)3) sm{63) + 2piyP2zE3 sin((/)3) sin(6'3) + 
2cos((/.3)2^2Pi.^3sin(e3)^ - 2cos(</.3)2£;ip2^^3sin(03)^ + 

(A.ll) 



2E2P1ZE3 sin(</.3)2 sin(03)^ - 2Eip2zE3 sin(</)3)2 sin(03)^ 



with X = 2p3.(pi + P2)/E3 and s standing for the squared invariant mass of the cohiding 
hadrons. 

If we treat the variables Ei and a = 2pi.p3 as two independent parameters, the 
components of the three- momentum Pi of the missing particle and the energy E3 = \p^\ of 
the massless visible particle can be expressed as the solution of the following linear system 
of four equations 

{Pi+P2)^ = si2 (A.12a) 

{Pl + P2 + PS)^ = Sl23 (A.12b) 

pix + E3 sin 63 cos </.3 = -p'^r'"'' (A. 12c) 

Ply + E3 sin 63 sin ^3 = -P^'^''''''' (A. 12d) 

that is parametrized by the momentum p2, by the angles ^3 and (p3, by the total transverse 
momentum ^branches q£ g^ij ^]^g branches represented, by the blobs in Figure 14 and by 



the variables a and Ei. The next step is to determine the values of the variables a and 
El. The mass-shell condition for the missing particle of momentum pi and the equation 
2pi-P3 = ct defines a system of two coupled quadratic equations in the variables Ei and 
a, parametrized by the momenta of the blocks. This system can be solved analytically. 
There are up to four solutions for Ei and a. Each solution that is physical {i.e., such that 
\p3\ > 0, El > and each of the Bjorken fractions qi,q2 is between and 1) corresponds to 
a distinct phase-space point at which the Jacobian in Eq. ( |A.ll ) and the integrand must 
be evaluated. 

A.4 MB D 

The notation for the phase-space variables associated with this constrained sector is given 



in Figure 15. The momenta of the missing particles are denoted by pi and p2, the momenta 
of the branches connected to the main block are denoted by p3, p4,, p^ and pQ. The Bjorken 
fractions are denoted by qi and 52- The variables Sij and Sjjfc refer to the invariants {pi+Pj)"^ 
and (pi+Pj+Pk)'^, respectively. The standard phase-space parametrization associated with 
this main block reads 
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Figure 15: Notation for the kinematics of MB D. 



The four-vectors Pjn and Pgn refer to the total momenta in the initial and final states, 
respectively. In our procedure we apply a change of variables that leads to the following 
parametrization of the phase-space measure 



The Jacobian J of this transformation is given by 

E1E2 



dsi3dSi34,dS25dS256 X J- 



J 



8s 



E3IE5 [p34z{PlyP2zP5ex " PlxP2zP56y 



-PlyP2xPmz +PlxP2yP56z) + Plz{-P2zP34yP56x + 
P2zP34xP56y - P2yP34xP56z + P2xP3AyP56z)] + 
{E^6P2z - E2P56z){PlzP34yP5x " PlyP34zP5x " PlzP34xP5y + 
PlxP3AzP5y) + [E56{PlzP2yP34x - PlzP2xP3Ay + PlyP2xP34z " 
PlxP2yP34z) + E2{pizP3AyP56x " PlyP34zP56x " PlzP34xP56y + 
PlxP3AzP56y)]P5zj + EsAE5P2z(pizP3yP56x - PlyP3zP5Gx 
-PlzP3xP56y +PlxP3zP56y) + E^{pizP2yP3x - PlzP2xP3y 
+PlyP2xP3z - PlxP2yP3z)P56z - {E^&P2z - E2P56z) 
iPlzP3yP5x - PlyP3zP5x - PlzP3xP5y + PlxP3zP5y) 
-[E56{PlzP2yP3x - PlzP2xP3y + PlyP2xP3z - PlxP2yP3z) + 
E2ipizP3yP56x - PlyP3zP56x - PlzP3xP56y + PlxP3zP56y)]p5z> 

El\[E5{p2z{-p3AzP3yP56x + P3AyP3zP56x + 
P3AzP3xP56y - P3AxP3zP56y) + 

{-P2yP3AzP3x + P2xP3AzP3y + P2yP3AxP3z - P2xP3AyP3z)P56z] + 
[E56P2Z - E2P5(iz){P3AzP3yP5x - P3AyP3zP5x - PSAzPSxPby + 
P3AxP3zP5y) + {E^6{P2yP3AzP3x - P2xP3AzP3y - P2yP3AxP3z + 
P2xP3AyP3z) + E2{p3AzP3yPmx - P3AyP3zP56x - P3AzP3xP56y + 



+ 



P3AxP3zP56y 



)]P5z} 



(A.14) 



(A.15) 
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where Eij = Ei + Ej , pij = pi + pj , and s is the squared invariant mass of the colhding 
hadrons in their center-of-mass frame. If we treat the variables Ei and £'2 as independent 
parameters, then the components of the three- momenta Pi,P2 of the missing particles can 
be expressed as the solution of the following linear system of six equations 



{pi + psY' 



Sl3 



(Pl +P3+P4.f = S134 

{P2+P5f = S25 
{P2+P5+P6f 



•S256 



^ _L ,r, — ^branches 
Plx + P2x — -PTx 



Ply + P2y 



^branches 
-PTy 



(A.16a) 
(A.16b) 
(A.16c) 
(A.16d) 
(A.16e) 
(A.16f) 



that is parametrized by the momenta p^, . . . ,pQ of the branches connected to the main 
block, by the total transverse momentum p™-anches q£ ^y[ the branches represented by the 



blobs in Figure 15 and by the variables Ei and £'2- The next step is to determine the 
values of the variables Ei and E2 ■ The mass-shell conditions for the two missing particles 
of momentum pi and p2 define a system of two coupled quadratic equations in the variables 
El and E2, that can be solved analytically. There are up to four solutions for Ei and £2- 
Each solution that is physical (i.e., such that E2 > 0, Ei > and each of the Bjorken 
fractions qi,q2 is between and 1) corresponds to a distinct phase-space point at which 
the Jacobian in Eq. ( |A.15 ) and the integrand must be evaluated. 



A.5 MB E 

The notation for the phase-space variables associated with this main blob is given in Fig- 



ure 16. The momenta of the missing particles are denoted by pi and p2, the momenta 





Figure 16: Notation for the kinematics of MB E. 

of the branches directly connected to these missing particles are denoted by p^ and p4. 
The Bjorken fractions are denoted by gi and q2- The variables Sij refer to the invariants 
{pi +Pj)'^ and s denotes the squared invariant mass of the colliding partons. The standard 
phase-space parametrization associated with this MB reads 



dqidq2 



cP'pi (fp2 

(27r)32£i (27r)32£2 



(27r)^5^ (Pin - Pi 



fin J 



(A.17) 
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The four-vectors Pin and Pgn refer to the total momenta in the initial and final states, 
respectively. In our procedure we apply a change of variables that leads to the following 
parametrization of the phase-space measure 

-dydsdsi^ds2i x J, (A. 18) 



IQtt'^EiE2 

where y is the rapidity of the colliding partons in the lab frame. The Jacobian J of this 
transformation is given by 

_ E1E2 
" 4s 

PlxP2yP3z) + E2PlzP3yPix " EiP2zP3yPix " E2PlyP3zP4.x + EiP2yP3zPix 

-E2PlzP3xP4y + EiP2zP3xP4y + E2PlxP3zPAy " EiP2xP3zP4y + (-E'2Plj/P3x + 

-EiP2yP3x - E2PlxP3y + EiP2xP3y)Piz + E2,{-pizP2yPix + PlyPlzPix 



J 



EA{pizP2yP3x - PlyP2zP3x " PlzP2xP3y + PlxP2zP3y + PlyP2xP3z 



+PlzP2xP4y - PlxP2zPAy - PlyP2xPAz + PlxP2yPAz 



(A.19) 



If we treat the variables £'1, E2 and p2y as independent parameters, the other compo- 
nents of the momenta pi , p2 of the missing particles can be expressed as the solution of the 
following linear system of five equations 

{pi+P3? = si^ (A.20a) 

{P2+PA? = S2A (A.20b) 

Pix+P2x = -p''J''''''"'' (A.20c) 

Ply+P2j; = -P;;'^"°'^''^ (A.20d) 

Plz + P2z = sinh(y)sl/2 _ ^branches (^206) 

that is parametrized by the momenta of the branches p^ and p4, by the total momentum 
^branches q£ g^^j ^j^^ branches represented by the blobs in Figure |l^, by the rapidity y and the 
invariant mass s^" of the colliding partons and by the variables -Ei, E2 and p2y The next 
step is to fix the values of the variables Si, E2 and p2y The variable Ei can be expressed 
as a linear function of E2 : 

El = COsh(y)sl/2 _ ^branches _ ^^^ ^^21) 

Then the mass-shell conditions for the two missing particles define a system of two coupled 
quadratic equations in the variables E2 and p2y In this case, the quartic terms of the two 
equations have the same coefficients, and the system reduces to a linear equation and a 
quadratic equation. There are up to two solutions for E2 and p2y Each solution that is 
physical (i.e., such that E2 > 0, Ei > 0) corresponds to a distinct phase-space point at 
which the Jacobian in Eq. ( [A.19 ) and the integrand must be evaluated. 



A.6 MB F 

The notation for the phase-space variables associated with this constrained sector is given in 



Figure 17, The momenta of the missing particles are denoted by pi and p2, the momenta 
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Figure 17: Notation for the kinematics of MB F. 



of the branches directly connected to these missing particles are denoted by p^ and ^4. 
The Bjorken fractions are denoted by gi and (72- The variables Sjj refer to the invariants 
{pi+Pj)"^- The standard phase-space parametrization associated with this main block reads 

d?pi d?P2 



dqidq2 



(2^)^5^ (Pi, - Pi 



finj 



(A.22) 



' (27r)32£;i {2tt)^2E2 

The four-vectors Pjn and Pa^ refer to the total momenta in the initial and final states, 
respectively. In our procedure we apply a change of variables that leads to the following 
parametrization of the phase-space measure 

1 



Wir'^EiEi 
The Jacobian J of this transformation is given by 

E1E2 



dqidq2dsi3ds24 x J- 



(A.23) 



J 



E4,{pizP2yP3x -piyP2zP3x - PlzP2xP3y + PlxP2zP3y + PlyP2xP3z- 

PlxP2yP3z) + E2PlzP3yPAx " EiP2zP3yP4.x " E2PlyP3zPAx + EiP2yP3zP4x 

E2PlzP3xPAy + EiP2zP3xP4y + E2PlxP3zPAy " EiP2xP3zP4y + {E2PlyP3x + 

-Eip2yP3x - E2PlxP3y + EiP2xP3y)PAz + E3{-pizP2yP4x + PlyP2zP4x 

-1 

+PlzP2xP4y - PlxP2zP4y " PlyP2xPAz + PlxP2yP4z 



(A.24) 



If we treat the variables Ei, E2 and p2y as independent parameters, the other compo- 
nents of the momenta pi , p2 of the missing particles can be expressed as the solution of the 
following linear system of five equations. 



(pi +P3f = Si3 
{P2 +P4f = S24 



„ !^ „ ^branches 

Plx + P2x — -Px 



(A.25a) 
(A.25b) 
(A. 25c) 
(A.25d) 

^^'^{qi - ^2)/2 - ^branches (^ 25e) 

that is parametrized by the momenta ps and p4, by the total momentum jf^^'^^^^^ of all 



Ply + P2y 
Plz + P2z 



branches 

y 



the branches represented by the blobs in Figure 17, by the Bjorken fractions gi, 52 and by 
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the variables Ei, E2 and p2y The next step is to fix the values of the variables i?i, E2 and 
P2y The variable Ei can be expressed as a linear function of E2: 

El = S^l'^{qi + q2)/2 - ^branches _ ^^_ ^^26) 



The mass-shell conditions for the two missing particles with momenta pi and p2 define a 
system of two coupled quadratic equations in the variables E2 and p2y In this case, the 
quartic terms of the two equations have the same coefficients, and the system reduces to 
a linear equation and a quadratic equation. There are up to two solutions for E2 and 
P2y Each solution that is physical (i.e., such that £'2 > 0, E'l > 0) corresponds to a 



distinct phase-space point at which the Jacobian in Eq. ( A. 24 ) and the integrand must be 
evaluated. 



B. Phase-space measure associated with the secondary blocks 

B.l SB A 

The notation for the phase-space variables associated with this secondary block is given in 
Figure ll^. The momentum of the missing particle is denoted by pi, the momenta of the 




Figure 18: Notation for the kinematics of SB A. 

three branches connected to the block are denoted by p2, Ps and p^. The variables S12, 
S123 and S1234 refer to the invariants (pi +^2)^, {pi + P2 +^3)^, and (pi + P2 + Ps +Pa)'^, 
respectively. The standard phase-space parametrization associated with this block reads 

(2^0^2^- ^^-^^ 

In our procedure we apply a change of variables that leads to the following parametrization 
of the phase-space measure 



1 



dsi2dsi23dsi234 X J. 



{2'Kf2Ei 

The Jacobian J of this transformation is given by 
El 

J = -^ EA{PlzP2yP3x -PlyP2zP3x - PlzP2xP3y + PlxP2zP3y + PlyP2xP3z 
o 

-PlxP2yP3z) + E2PlzP3yP4x " EiP2zP3yPAx - E2PlyP3zPAx + EiP2yP3zPAx 

-E2PlzP3xP4y + EiP2zP3xP4y + E2PlxP3zP4y " EiP2xP3zP4y + {E2PlyP3x 

-EiP2yP3x - E2PlxP3y + EiP2xP3y)PAz + E3{-pi-^p2yPAx + PlyP2zPAx 

+PlzP2xPAy - PlxP2zPAy " PlyP2xPAz + PlxP2yP4z 



(B.2) 



(B.3) 
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If we treat the variable Ei as an independent parameter, the components of the three- 
momentum Pi of the missing particle can be expressed as the solution of the following 
linear system of three equations 

{pi+P2f = si2 (B.4a) 

{Pl + P2 + Psf = Sl23 (B.4b) 

(Pl +P2+P3+ Pif = •S1234 (B.4c) 

that is parametrized by the momenta p2 iPSi Pa of the branches connected to the block, 
and by the variable Ei . The next step is to fix the value of the variable Ei . The mass-shell 
condition for the missing particle with momentum pi defines a quadratic equation in the 
variable Ei. There are up to two solutions for Ei. Each solution that is physical (i.e., 
such that El > 0) corresponds to a distinct phase-space point at which the Jacobian in 
Eq. (|B.3|) and the integrand must be evaluated. 



B.2 SB B 

The notation for the phase-space variables associated with this secondary block is given in 
Figure ^. The momentum of the missing particle is denoted by pi, the momenta of the 



9 

/P3 



S123 S12 




Figure 19: Notation for the kinematics of SB B. 

two branches connected to the block are denoted by p2 and p^. The variables si2 and S123 
refer to the invariants (^1+^2)^! and (^1+^2+^3)^, respectively. The standard phase-space 
parametrization associated with this block reads 

^'^^ (B.5) 



(27r)32£;i ■ 

In our procedure we apply a change of variables that leads to the following parametrization 
of the phase-space measure 

1 



(27r)32£;i 



!)i(isi2(isi23 X J, (B.6) 



where (pi denotes the azimuthal angle of particle i. The Jacobian J of this transformation 
is given by 



T ^1 



COs((/)i - (l)2)E3P2TPlz + COs((/)i - (j)^) E2P3TPIZ + E3P1TP2Z 

-1 
COs{(j)l - (l)3)EiP3TP2z - E2PITP3Z + COs{(pl - (l)2)EiP2TP3z ■ (B.7) 
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If we treat the variable Ei as an independent parameter, the transverse momentum piT 
and the momentum component piz of the missing particle can be expressed as the solution 
of the following linear system of two equations 

{Pi+P2f = si2 (B.8a) 

{Pl + P2 + P3f = Sl23 (B.8b) 

that is parametrized by the momenta p2 and pa of the branches connected to the block, 
by the azimuthal angle 0i and by the variable Ei . The next step is to fix the value of the 
variable Ei . The mass-shell condition for the missing particle with momentum pi defines a 
quadratic equation in the variable Ei . There are up to two solutions for Ei . Each solution 
that is physical (i.e., such that Ei > 0) corresponds to a distinct phase-space point at 
which the Jacobian in Eq. (|B.7D and the integrand must be evaluated. 



B.3 SB C/D 

The notation for the phase-space variables associated with this secondary block is given in 
Figure pQ. The momentum of the missing particle is denoted by pi , the momentum of the 




Sl2 

bil 

Figure 20: Notation for the kinematics of SB C/D. 

branch connected to the block is denoted by p2- The variable si2 refers to the invariant 
{pi +^2)^- The standard phase-space parametrization associated with this block reads 

^'^^ (B.9) 



(27r)32£;i ■ 



In our procedure we apply a change of variables that leads to the following parametrization 
of the phase-space measure 

^ ' hd0idsi2 X J, (B.IO) 



{2tt)^2Ei 

where 61 and (pi denote the polar and azimuthal angles of the missing particle. The 
Jacobian J of this transformation is given by 

1 



J = — sme^ilpil 



\P1\E2 - Eip-^.P2 



(B.ll) 



If we treat the variable Ei as an independent parameter, the momentum modulus \pi\ of 
the missing particle can be expressed as the solution of the following linear equation 

ipi+P2f=Sl2 (B.12) 



37 



that is parametrized by the momentum p2 of the branch connected to the block, by the 
polar and azimuthal angles 9i, (pi, and by the variable Ei. The next step is to fix the 
value of the variable Ei . The mass-shell condition for the missing particle with momentum 
pi defines a quadratic equation in the variable Ei. There are up to two solutions for Ei. 
Each solution that is physical {i.e., such that Ei > 0) corresponds to a distinct phase-space 
point at which the Jacobian in Eq. (B.ll) and the integrand must be evaluated. 



B.4 SB E 

The notation for the phase-space variables associated with this secondary block is given in 
Figure ^. The momenta of the visible particles are denoted by pi and p2, the momentum 




S123 



Sl2 




Pi 



P2 



Figure 21: Notation of the kinematics for SB E. 



of the branch connected to the block is denoted by p^. The variables S12 and S123 refer to 
the invariants (pi +^2)^, and (pi + P2 +^3)^- The standard phase-space parametrization 
associated with this block reads 



d^pi (fp2 



(B.13) 



(27r)32Ei (27r)32^2 ' 

In our procedure we apply a change of variables that leads to the following parametrization 
of the phase-space measure 

1 



-d9id(j)id92d(f>2dsi2dsi23 x J, 



(B.14) 



{27r)^4:EiE2 
denote the polar and azimuthal angles of particle i. The Jacobian J of 



where 9i and 

this transformation is given by 

J = —r\Pi\ sin0isin02 



\Pi\\P2\/Ei - \P2\fi2){E3 - \P3\f23) 



-{E^Ml/Ei - \p^\fi3){Ei - fi2\Pi\ 



(B.15) 



where fij stands for Pi-Pj/\Pi\\Pj\- The values for the momenta \pi\ and IP2I can be 
obtained by solving the following linear system of equations 



{Pl+P2f 
{Pl +P2+P3f 



S12, 

Sl23- 



(B.16a) 
(B.16b) 



By subtracting Eq. ( |B.16b ) from Eq. ( |B.16a ), we obtain an expression for Ei that is a first 
order polynomial in \pi\ and \P2\- Inserting this expression into Eq. ( [B.16a| ) and into the 
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equation defining the mass-shell condition for the particle of momentum pi, we obtain a 
system of two quadratic equations in |p;^| and IP2I parametrized by the momentum ^3, by 
the invariants S12, S123 and by the angles 9i, 62, (pi and 4>2- This system can be solved 
analytically. There are up to four solutions for the modulus |p^| and \P2\- Each solution 
that is physical (i.e., such that \pi\ > and IP2I > 0) corresponds to a distinct phase-space 
point at which the Jacobian in Eq. ( B.15| ) and the integrand must be evaluated. 
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• ' Abstract: Matrix element reweighting is a powerful experimental technique widely em- 

f^ ' ployed to maximize the amount of information that can be extracted from a collider data 

^P ■ set. We present a procedure that allows to automatically evaluate the weights for any 

process of interest in the standard model and beyond. Given the initial, intermediate and 
final state particles, and the transfer functions for the final physics objects, such as leptons, 
jets, missing transverse energy, our algorithm creates a phase-space mapping designed to 
C^ \ efficiently perform the integration of the squared matrix element and the transfer functions. 

The implementation builds up on MadGraph, it is completely automatized and publicly 
available. A few sample applications are presented that show the capabilities of the code 
and illustrate the possibilities for new studies that such an approach opens up. 
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1. Introduction 

Along with the ongoing experimental activity in Run II at the Tevatron and at the now 
operational Large Hadron Collider at CERN, in the last years a significant effort by the 
high-energy community (including both theorists and experimentalists) has been devoted 
to devise new and more efficient strategies to identify physics beyond the standard model 
in collider data. In many of the new physics scenarios, new states exist at TeV scale that 
decay very quickly, and are not expected to leave any trace that can be reconstructed in 
the detector. Hence their existence and properties must be inferred from the distributions 
and properties of standard model particles that can originate from the decay of heavier 
not-yet-discovered resonances. 

The problem of identifying such decay patterns and from those of measuring the prop- 
erties of the new states is particularly intricate when the expected experimental signatures 
involve a complex final state, typically with several jets, leptons and missing energy. The 
latter, in particular, characterizes many of the models that aim at providing a candidate 
for dark matter consistent with the present observations. Several methods have been de- 
veloped during the last few years to improve on and eventually overcome this difficulty. 
For the purpose of identifying new physics, it is quite natural to first consider an approach 
that is not biased by strong theoretical assumptions, as the current knowledge of the vi- 
able theories that could lead to the production of new particles, is somehow limited. In 
this context, different methods have been proposed to measure the mass spectrum of the 
new states in a model-independent way: specific observables are suggested and built that 
are mostly sensitive to the masses of the new heavy resonances entering the decay chains. 
The final power of a given method is a balance between how well the information of the 
visible quantities is exploited to constrain the unknown masses and the dependence on the 
experimental and theoretical systematic uncertainties. Examples in the literature include 
the end-point method — that is based on the end-point regions of the invariant mass dis- 
tributions built from visible particles — , and the polynomial method — that attempts to 
reconstruct the whole event from the visible momenta — . These two methods can also be 
combined to give a better constraint on the mass spectrum. The Mt2 method [?, ?], and 
its generalizations/recent developments — though based on more complicated observables- 
also follow the same philosophy. For a recent review on these kinematics methods see [?] 
and references therein. 

These model-independent methods will be determinant in constraining the mass spec- 
trum of new resonances. However, hy construction, most of them will not exploit or provide 
any information on other properties of the new particles, such as spin and coupling struc- 
ture. As another example, the precise measurement of the absolute mass of each particle 
entering a specific decay chain that ends with two missing particles remains challenging, 
especially in the case of short-length decay chains. In this context, it is useful to consider 
complementary and model-dependent tools for the investigation of properties of the new 
physics states. 

Matrix element reweighting [?, ?, ?, ?, ?, ?, ?] is an example of such a method. It 
dramatically differs from the previous ones in that it requires at least one theoretical as- 



sumption as a starting point. Each assumption specifies the rules needed to compute the 
probability distribution associated with the process under study. Following a Bayesian ap- 
proach, the method assigns a probability to each hypothesis given a sample of experimental 
events, and in this way provides a discriminator among the different hypotheses. Another 
important feature of the method is that it makes maximal use of both experimental in- 
formation and the theoretical model (via the amplitude) on an event-by-event basis. This 
optimal use of the experimental data as well as the theoretical knowledge opens the door 
to new studies not only on mass measurements but also on the identification of the spin 
and coupling type of new particles. 

The matrix element method has been extensively exploited in the last decade. The 
best known example is its application to top-quark-pair production — investigated both by 
the CDF and D0 collaborations [?, ?, ?, ?, ?] — which has led to the single most precise 
measurement of the top-quark mass. Recently, it has contributed to the observation of 
single-top production [?, ?] and to set an upper limit on the boson production of a standard 
model Higgs [?], currently excluded in the mass region 158 GeV < mn < 175 GeV at 95% 
confidence level [?]. 

In principle, the matrix element method is expected to provide a powerful investigation 
tool in many other analyses, in particular those dedicated to the search of new resonances 
and the study of their properties. However, in practice, its application is not straightfor- 
ward. In order to evaluate the weights to be attached to each experimental event, a difficult 
convolution of the theoretical information on the hard scattering (i.e., the matrix element 
squared) with the experimentally available information on the final state (encoded in the 
so-called transfer functions) has to be undertaken. The numerical efficiency (and therefore 
the speed) of such integration is currently a serious limitation. The matrix element squared 
as well as the transfer functions present variations by several orders of magnitude in differ- 
ent regions of the phase space. To overcome this difficulty, the integration technique has 
to be efficiently adapted to the shape of the integrand. To our knowledge, this problem 
has only been solved in very specific cases. 

In this work, we propose a general algorithm aimed at evaluating the weights appearing 
in the matrix element method. Given an arbitrary decay chain and the associated transfer 
function, our procedure first automatically assigns the optimized phase-space mappings 
designed to match as much as possible the peaks in the integrand, and then performs the 
phase-space integrations to evaluate the weights. Our implementation, which is fully auto- 
matic, is based on MadGraph, as it uses its matrix element amplitudes and the information 
on the topology of the diagrams. We dub the corresponding public code MadWeight. 

The paper is organized as follows. In Section |2| we review the basic features of the 
matrix element method. In Section ^ we expose our algorithm for the computation of the 
weights in the matrix element method. We present some illustrations in Section ^ and our 
conclusion in the last Section. 

2. The matrix element method 

As mentioned in the introduction, the matrix element method is a procedure to extract 



theoretical information, in the form of set of parameters a from a sample of experimen- 
tal events.^ Let us identify an event by the set x of experimentally available quantities 
(such as transverse momenta, rapidities, and so on). For each observed event a conditional 
probability P{x\Q.),i.e., a weight, is built that quantifies the "agreement" between the the- 
oretical framework a and the experimental event x. In the computation of the weights, one 
factorizes high-energy effects associated with the production of a parton-level configuration 
y into a calculable probability Pa{y). The evolution of the parton-level configuration y 
into a reconstructed event x in the detector is modeled by a transfer function W{x, y). As 
a result, the weight of a specific event x is of the form 

P{x\cx) = j dyP^{y)W{x,y). (2.1) 

In the specific case of a hadron collider, the parton-level probability Pa{y) can be expressed 
as a product of the squared matrix element |MQ,p(y), the parton distribution functions 
(pdf's) fi{qi) and f2{<l2) and the phase-space measure d^{y), such that the weight reads 

P(a,|a) = 1- [ d^{y)dqidq2h{qi)f2{q2)\M^\\y)W{x,y) . (2.2) 

0"a J 



The normalization by the total cross section a a in Eq. ( |2.2| ) ensures that P{x\oi.) is a 
probability density^: j P{x\cx)dx = 1. Once this probability density has been computed 
for each event Xi, the most probable value for a can be obtained through a likelihood 
maximization method. 

Eq. ( p.2| ) is central to this paper as it provides an explicit definition of the weight to be 
associated with a given event in terms of the convolution of tree- level matrix element, the 
pdf's and the transfer functions. One of the main working assumptions in the application 
of the matrix element method is that the transfer functions are " factorisable" ,i.e., they can 
be written as the product of single-particle resolution functions 

n 

W{x,y) = J{Wi{x\y'), (2.3) 

1=1 

where x* and y* stand for the measured quantities and the phase-space variables associated 
with the particle i, respectively. In practice, a further simplification is employed, where the 
transfer function associated with a single reconstructed object (such as a jet or a lepton) 
is written as a product of resolutions associated with the physical quantities measured in 
the detector: 

W,{x\y') = W^{x\y')Wl>{x\y')Wt{x\y'), (2.4) 



^Normally a. labels the different parameters in a given model (such as, for example, a mass or the value 
of a coupling). In this paper, however, we use a more general definition that also includes labelling different 
physics models. 

^We assume that the transfer function is also normalized to 1. 



where E, rj and (j) ^-re the reconstructed energy, rapidity and azimuthal angle, in most of 
the general purpose detectors the direction of a visible particle^ is well measured, so that 
the associated transfer function can be modeled by a narrow Gaussian. On the other hand, 
the resolution in energy strongly depends on the particle's type. For leptons it can be 
taken as a narrow Gaussian function whereas for jets a more involved parametrization of 
the resolution function is needed. 

In this work we provide a general solution to the problem of performing the integration 
in Eq. ( |2.2| ) in an efficient way. To better grasp the challenge that computing the integral 
in the numerator of Eq. ( |2.2| ) poses, it is useful to consider two limiting cases, where the 
problem simplifies. 

First let us imagine to have an "ideal" detector that could measure exactly the energies 
and momenta of all final state particles (including normally invisible ones), i.e., W{x, y) = 
6{x — t/).^ In this case no integration would be necessary and the weight in the numerator 
would be proportional to the corresponding squared matrix element, \Ma\'^{x). Nowadays, 
the determination of |Mo,p(a;) at the tree-level can be done automatically by several public 
codes and poses no difficulty. So apart from the normalization, discussed below, the weight 
calculation would therefore be trivial. 

As a second limiting case, one can also consider an ideal "no detector" option, i.e., 
choose the transfer function W{x,y) = 1. Then the integration would reduce to the 
computation of the total cross section, i.e., that of the denominator of Eq. (|2.2| ) as P{x\a) = 
1. This problem is not an easy one on its own: the matrix element has a very complicated 
peak structure, corresponding to the propagators of the Feynman diagrams being large. 
However, by observing that the leading peaks come from the sum of the squares of each 
diagram, together with the fact that it is always possible to find a parametrization of the 
phase space in terms of invariants that maps exactly those in the propagators [?] , makes the 
problem treatable (see for example Refs. [?, ?] and the discussion in the following Section). 

For a realistic detector the situation is in between the two above, where some par- 
ticles are well measured (charged leptons), other less (jets), and some completely missed 
(neutrinos). In this case the integration becomes extremely difficult as it involves an inte- 
grand with simultaneous peaks in sets of different variables that it is not possible, even in 
principle, to disentangle. 

3. Computation of the weights 

The evaluation of multi-dimensional integrals is often approached by standard adaptive 
Monte Carlo techniques. These techniques are well illustrated in the computation of total 
cross sections: phase-space mappings that "flatten" specific peaks in the integrand are 
combined together in a multichannel integration. Here we also follow this approach. In our 



'^Throughout the paper, the expression visible particle refers to a lepton or a jet of which momentum is 
reconstructed in the detector. 

*For jets, the transfer functions also include genuine QCD effects like showering and hadronization. 
For the sake of the argument, we consider also "ideal" jets where the identification jet/parton is perfectly 
unambiguous. 



case, the phase-space mappings optimized for the computation of the weights in Eq. (^ 
are rather involved because of the complex structure of peaks in the integrand, as it was 
discussed in the previous section. 

In this section we present our integration procedure and its implementation in a fully 
general algorithm. We first recall the basic principle of an adaptive Monte Carlo integration 



in Section 3.1 and then we describe the phase-space mappings optimized for the computa- 



tion of the weights in Section 3.2. We explain how we build a phase-space generator based 



on these new phase-space mappings in Section B^ and how we combine different phase- 



space mappings in a multi-channel integration in Section |3.4 We validate our phase-space 
generator with several checks in Section |3.5|. 



3.1 Adaptive Monte Carlo techniques 

Adaptive Monte Carlo integration is a powerful numerical technique for the integration of 
a highly non-uniform function. It consists of sampling randomly the volume of integration 
according to a probability density that is adjusted iteratively to the shape of the integrand. 
The probability density is parametrized by a separable function 

p{z)=pi{z^)p2{z^) ...Pdiz'') (3.1) 

where each factor pi is a step function. If such a parametrization of the probability density 
function is appropriate to approximate the shape of the integrand, the adaptive integration 
procedure speeds up the convergence by increasing the density of evaluations in the regions 
where the integrand is large. In the case of a very sharp integrand, this condition is 
essentially fulfilled provided that the strength of each narrow peak in the integrand is 
associated with a single variable that in turn can be mapped onto one variable of integration 
z*. In that case, the integrand expressed in the parametrization z is of the form 



fi^)= Uhi^')] ^Riz) (3.2) 




where the functions /j's may vary abruptly while the "remainder" non-factorisable function 
R{z) is essentially flat over the region under integration. 

If the integrand expressed in the phase-space mapping z presents a structure of sharp 
peaks that does not follow the factorized form in Eq. (^.2|), the adaptive integration pro- 
cedure is bound to fail. However, if enough information about the shape of the integrand 
is available, a first change of variables z ^ z' = P{z) that rotates the axes of integration 
can sometimes be applied such that in the new phase-space mapping z' , the importance 
of each peak in the integrand is controlled by a single variable of integration. After this 
change of variables is applied, the integrand expressed in the new variables z' is of the form 
given by Eq. (|3.2| ), and the separable density function p{z') can be successfully adapted to 
the shape of the integrand. 

We will use the adaptive Monte Carlo integrator VEGAS [?] to carry out the inte- 
gration in Eq. (p. 21). Thus the efficiency in computing the weights will depend on the 



parametrization of the phase-space measure that is used in the adaptive Monte-Carlo in- 
tegration. The optimized phase-space mappings are such that for each narrow peak either 
in the transfer function or in the matrix element, the variable that controls the strength 
of that peak is mapped onto a single variable of integration in the parametrization of the 
phase-space measure, in which case the integrand expressed in that parametrization has 
the form given in Eq. (^Z 



3.2 The ne^v phase-space mappings 

For the computation of the weights, there is generally no simple phase-space parametriza- 
tion that maps all the peaks in the integrand and in which the boundaries of the phase- 
space volume can be easily expressed. Our strategy is to start from the following standard 
parametrization of the phase-space measure 

^,is [tJ \Pi\'^d\Pi\ sin 9ideid(l)i\ /o ^4x4 / , V^ I a o\ 

'^^=111 2E(2tt)^ j dqidq2i27r) 6 \pi+P2-2^Pj\, (3.3) 

where i = 3, . . .n labels the final particles. In this parametrization, the strength of each 
peak in the transfer function is already mapped onto a single variable of integration, whereas 
none of the propagator enhancement in the squared amplitude is. Identifying the Lorentz 
invariants associated with the Breit-Wigner resonances and expressing them as functions of 



the integration variables in Eq. (3^) is straightforward. The difficult task is then to invert 
these functions in order to derive a phase-space measure that is parametrized by both 
these Lorentz invariants and the variables mapping the peaks in the transfer function. 
Along with this inversion, the 6 function associated with energy-momentum conservation 



in Eq. (3^) has to be integrated out. The resulting phase-space mappings can then be 
used in an adaptive Monte Carlo integration to compute the weights. 

These optimized phase-space mappings can be defined by specifying the transformation 
of the phase-space measure parametrization in Eq. ( |3.3| ) from which they result. So in this 
Section, we will describe the expression of this transformation in a generic case, as it is a 
convenient way to introduce the new phase-space mappings. For an arbitrary process, the 
transformation that leads to the appropriate parametrization of the phase-space measure 
can be carried out by organizing the integration variables in the standard parametrization 



in Eq. (3^) into different subsets of variables to which a suitable change of variables is 
applied. Each subset of variables and its associated change of variables will be called a 
block in the following. 

The first phase-space block that needs to be identified is called the main block (MB), 
and it includes some of the integration variables appearing in Eq. ( |3.3| ) to which a trans- 
formation is applied so that the 6 function associated with energy-momentum conservation 
is integrated out. The same transformation may also map some invariants entering in the 
expression of specific propagators to new variables of integration in the expression of the 
phase-space measure. The identification of the main block and the form of the associated 
transformation of variables is discussed in the following Section. The integration variables 
appearing in Eq. (|3.3|) that do not belong to the main block also experience a transforma- 



tion that can be expressed in terms of secondary blocks, as explained in Section 3.2.2 
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Figure 1: Illustration of a decay chain with no missing particle: (a) the full topology, (b) the 
corresponding reduced diagram. The variables in the main block are written explicitly. 



3.2.1 Identification of the main block 

The main block includes a certain number of integration variables among those appearing 
in Eq. ( [3.31 ). These variables are adjusted as a function of all other kinematic quantities 
associated with the decay chain to enforce the conservation of energy and momentum. 

We start by discussing the choice of the main block in the case of two specific decay 
chains, and then generalize to the case of an arbitrary decay chain. We first consider 
a topology with no missing particle. An example of such a decay chain is illustrated in 
Figure |l|(a) . In that case, it is natural to include the initial proton momentum fractions of 
partons (called Bjorken fractions) in the main block, as the integrand does not show any 
sharp sensitivity in these variables. The angle of any visible particle should be excluded 
from the main block as it controls the strength of a narrow resolution function W^ or W^ 



in Eq. (2.4): it should therefore be maintained as an integration variable to ensure the 
integration of the associated peak. However, normally the resolution in energy is much 
poorer than the resolution in angles. For this reason, a relatively efficient choice is to 
complete the main block by adding two momentum variables \pj\ and |p,| of particles i 
and j that are relatively less constrained by the transfer function. In this example, the 
main block contains exactly four variables, and the effect of the variable transformation is 
to integrate out the 5 function with these four variables. 

We next move on to the case of a decay chain including missing particles in the final 
state. The phase-space variables associated with the momenta of missing particles are not 
directly constrained by the transfer function. Therefore they do not need to be mapped 
onto variables of integration in the phase-space mapping. We can identify the main block 
by selecting the momentum components of some missing particles instead of the energies 
of visible particles. A specific example of topology with missing particles is displayed in 
Figure §(a). One way to define the main block is to choose the set including the Bjorken 
fractions and the momentum components of the missing particle shown as a thick line in 



Px,Py,Pz 



(a) (b) 

Figure 2: Illustration of a decay chain with two missing particles (identified by the letter v): (a) 
the full topology, (b) the corresponding reduced diagram. The initial (resp. final) variables of the 
transformation associated with the main block are written explicitly, and the corresponding legs 
are shown as thick lines (resp. dashed lines). 

Figure 0(a). The change of variables associated with this main block remove these five 
variables from the set of integration variables in order to integrate out the 6 function in 
Eq. ( |3.3D and to map the invariant mass of the resonance decaying into the missing particle 
onto a variable of integration in the new parametrization of the phase-space measure. 

In order to generalize the discussion of the choice of the main block to the case of an 
arbitrary decay chain, it is useful to introduce the following representation of the main 
block and the corresponding transformation of variables: 

• In a branch of legs with no kinematic variable in the main block, the decay products 
of the initial particle in the branch are shrunk into a blob. 

• The variables in the main block are written explicitly and the corresponding legs are 
shown as thick lines. 

• The new integration variables resulting from the change of variables associated with 
the main block are also written explicitly, and the corresponding intermediate legs 
are shown as dashed lines. 

• All other intermediate legs that do not touch a blob are hidden behind a rectangular 
box. 

We refer to the resulting graph as the reduced diagram. As an illustration, the reduced 
diagrams for the two topologies shown in Figures |l](a) and Q(a) are displayed in Fig- 
ures |l](b) and |2|(b), respectively. In general a blob in a reduced diagram may hide a 
complicated branch of particles. But the change of variables associated with the main 
block is parametrized only by the total momentum of each blob, it does not depend on the 
structure inside the blob. 

We use these reduced diagrams to represent the main block in general. The minimum 
number of variables in the main block is four. After the variable transformation associated 
with the main block is applied, the 6 function in Eq. ( ^.3| ) is integrated out with these 
four variables that therefore do not appear in the new phase-space mapping resulting 



from this transformation. The main block may contain p > 4 integration variables. The 
transformation that is applied in that case removes all these p variables from the set of 
integration variables appearing in the parametrization of the phase-space measure, and 
introduces p — 4 new variables of integration. Each of these new variables map a Lorentz 
invariant that controls the strength of a specific propagator in the matrix element. Thus the 
variable transformation associated with the main block not only enforces the conservation 
of total energy and momentum, but also may possibly optimize the parametrization of the 
phase-space measure for the integration of some specific Breit-Wigner enhancements. 

As the variables in the main block are removed from the phase-space mapping af- 
ter the corresponding transformation is applied, an integration variable in the standard 
parametrization that controls the strength of a narrow peak in the transfer function is 
preferentially not included in the main block, otherwise the phase-space mapping after 
transformation would loose track of this variable and would be inappropriate for the inte- 
gration of the corresponding peak. From this observation, it is clear that the choice of the 
main block will act upon the efficiency of the Monte Carlo integration. 

Each of the main blocks treated in our code is illustrated by a reduced diagram in 
Figure ^. Their number is restricted because we only keep the main blocks for which the 
corresponding change of variables is invertible analytically. The corresponding formulas 
are discussed in Appendix ^. 

MB A. The transformation removes the Bjorken fractions qi and q2 and the norm of the 
three-momenta Pj , Pj of two visible particles from the set of integration variables in 
the parametrization of the phase-space measure. 
Example: pp — )• ZZ — )• 4j. 

MB B. The transformation removes the Bjorken fractions qi and g2 and the 3-momentum 
of a missing particle from the set of integration variables in the parametrization of 
the phase-space measure. The new integration variable is the invariant mass of the 
particle decaying into the missing particle. 
Example: pp — t- Z{W~^ — )• l^u). 

MB C. The transformation removes the Bjorken fractions qi and q2, the 3-momentum 
of a missing particle and the energy of a massless visible particle^ from the set of 
integration variables in the parametrization of the phase-space measure. The new 
integration variables are the Lorentz invariants m* and m* associated with the 
mother particles decaying into the missing and the massless particles, respectively. 
Example: pp — )• [i — )• h{W'^ — )• l~^v)]\i ^ b{W~ — )• jj)] with massless b quarks. 

MB D. The transformation removes the Bjorken fractions qi and 52 and the 3- momenta 
of two missing particles from the set of integration variables in the parametrization 
of the phase-space measure. The new integration variables are the Lorentz invariants 
m* , m* , m* and m* associated with the first and second mother particles of each 



''In the case the particle is massive, the corresponding change of variables turns out to be not analytically 
invertible. 
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Figure 3: The reduced diagrams representing the six main blocks that have been investigated in 
our procedure. 



missing particle. 
Example: pp — )• [i — )• b{W^ 



l+u) 



b{W- 



l-iy)]. 



MB E. The transformation removes the 3-momenta of two missing particles from the 
set of integration variables in the parametrization of the phase-space measure. The 
new integration variables are the Lorentz invariants m* and m* associated with the 
mother particles of each missing particle. The integration over the Bjorken fractions 
is expressed as an integration over the invariant mass and the rapidity of the colliding 
partons. 
Example: pp ^ H ^ {W+ -^ l+i^){W- -^I'v). 

MB F. The transformation removes the 3-momenta of two missing particles from the 
set of integration variables in the parametrization of the phase-space measure. The 
new integration variables are the Lorentz invariants m^ and rrv^ associated with the 
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Figure 4: A example of reduced diagram for wlricli the transformation that is applied to the MB 
cannot be inverted by means of analytical formulas. 



mother particles of each missing particle. 
Example: pp — )• (W^ — )• l~^i')(W~ — )• Z~z/). 

For a given decay chain, the main block can be chosen in several different ways, possibly 
with some more efficient than others. Roughly speaking a specific main block is appropri- 
ate for the computation of the weights provided that it does not contain a variable that 
controls the strength of a very sharp resolution function, and provided that no very sharp 
Breit-Wigner distribution is included in the square of the reduced diagram. This second 
condition comes from the fact that none of the invariants entering into the expression of 
the propagators inside the box is mapped onto a single variable of integration in the new 
parametrization of the phase-space measure, so that parametrization is not appropriate for 
the integration of the corresponding propagator enhancements. 

It should be stressed once again that each transformation of variables that is applied to 
the variables in main block has been implemented in the code analytically: for an arbitrary 
phase-space point, given the momenta of all the legs ending by the blobs and the invariant 
mass of each leg represented by a dashed line in the reduced diagram, the variables in the 
main block are determined by means of analytical expressions (see appendix |^) . We have 
not explored further the possibility to use numerical procedure for this step. So any change 
of variables that is not invertible analytically has been excluded in our algorithm. 

For example, the main block displayed in Figure ^ with three missing particles has not 
been considered. In principle, the 3-momenta of the three missing particles and the Bjorken 
fractions could be adjusted to satisfy eleven constraints induced by the seven resonances 
and the conservation of 4-momentum. As this adjustment cannot be done by means of 
analytical expressions, this case is dealt with the main block B or C in our procedure. 

3.2.2 Identification of the secondary blocks 

Once the main block has been defined and the corresponding transformation applied, the 
parametrization of the phase-space measure associated with the m external legs in all the 
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Figure 5: The four secondary blocks with their corresponding change of variables. Initial (resp. 
final) variables are written explicitly, and the corresponding legs are represented by thick lines (resp. 
dashed lines). 



blobs of the reduced diagram is still the standard one: 



n 



|Pjpd|pj| svnOidOi 
2E,{2TTf 



(3.4) 



The parametrization of the phase-space measure is further transformed by organizing the 
integration variables in Eq ( ^.4| ) into secondary blocks, i.e., into subsets of variables, each 
of them being subject to a specific change of variables. The change of variables of the 
simplest block is just the identity, in which case the variables in this block are maintained 
in the parametrization of the phase-space measure in Eq. (f3.4|). For the purpose of listing 
the other changes of variables that we have investigated, it is useful to represent a block 
and its corresponding change of variables by a diagram in the following way: 

• The variables involved in the transformation are written explicitly. The legs asso- 
ciated with the initial variables appear as thick lines. The legs associated with the 
final variables -which correspond to the invariants that enter into the expression of 
specific propagators- are shown as dashed lines. 

• A blob stands for a branch of legs of which total momentum parametrizes the change 
of variables related to the block. 

In this representation, a blob can a priori be itself decomposed into several secondary 
blocks. However, as in the case of the main block, the change of variables associated 
with a given secondary block is only parametrized by the total momentum of each branch 
represented by a blob, it does not depend on the details of these branches. The number of 
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Figure 6: Illustration of the structure in blocks optimizing the parametrization of the phase-space 
measure in the case of a specific decay chain. The missing particles are indicated by the Greek 
letter i/. 

implemented blocks in our algorithm is reduced by our requirement of analytically invertible 
changes of variables. These blocks are displayed in Figure |5[ The changes of variables 
associated with each block are discussed in Appendix pi 

SB A. The transformation removes the 3-momentum of a missing particle from the set 
of integration variables in the parametrization of the phase-space measure. The new 
integration variables are the Lorentz invariants m* 
first, second and third mother particles of this missing particle. 



m* and m* associated with the 



SB B. The transformation removes the energy and the polar angle of a missing particle 
from the set of integration variables in the parametrization of the phase-space mea- 
sure. The new integration variables are the Lorentz invariants m* and m* associated 
with the first and second mother particles of this missing particle. 

SB C/D. The transformation removes the energy of a missing particle from the set of 
integration variables in the parametrization of the phase-space measure (version C). 
The new integration variable is the Lorentz invariant m*^ associated with the mother 
particle of this missing particle. In version D of this block, the missing particle is 
replaced by a visible particle, but the transformation remains the same one. 

SB E. The transformation removes the momenta \pi\ and \p2\ of two visible particles pro- 
duced by the same resonance from the set of integration variables in the parametriza- 
tion of the phase-space measure. The new integration variables are the Lorentz in- 
variants 7T1,*, and m,*„ associated with the first and second mother particles of these 
visible particles. The corresponding change of variables is invertible analytically only 
if at least one of the two visible particles is massless. 

This completes the description of the blocks that can be used in our procedure to op- 
timize the parametrization of the phase-space measure for the computation of the weights. 
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As a example, we illustrate a composition in blocks in Figure y in the case of a specific 
decay chain. 

3.3 One-channel phase-space generator 

Given an optimized phase-space mapping defined by its structure in blocks, one can then 
consider a phase-space generator built upon this phase-space mapping. The generation of 
an arbitrary phase-space point proceeds in two steps: 1) the generation of the integration 
variables appearing in the optimized parametrization of the phase-space measure, 2) the 
determination of the momentum of each leg in the decay chain and the computation of the 
jacobian factors. 

Concerning the first step, any variable of integration associated with the new phase- 
space mappings introduced in the previous section enters into one of the three following 
categories: 

1. The variable controls the strength of a resolution function. If the resolution function 
is a (5 distribution, the variable is fixed to the value associated with the experimental 
event. Otherwise, the grid of VEGAS is adapted such that the variable is generated 
according to a probability density that reproduces approximately the shape of the 
resolution function. 

2. The variable controls the strength of a propagator enhancement. In this case, the 
variable can be generated according to a probability density that reproduces exactly 
the shape of the propagator by using the inverse primitive function of a Breit-Wigner. 

3. The variable is either the polar or the azimuthal angle of a missing particle. In this 
case, the variable is generated according to a uniform distribution in the interval [0, vr] 
or [0, 27r] at the first iteration. The grid is adapted at each iteration to approximate 
the optimal probability density. 

Once the integration variables have been generated, the kinematics of the whole decay 
chain and the jacobian factors are computed. For each block, the formulas that give 
the expression of the external momenta as a function of the variables of integration are 
discussed in the Appendix. These formulas are parametrized by the momentum of the 
branches represented by the blobs that appear in the graphical representation in Figures ^ 
and ^. For this reason, one needs to fill the kinematic variables in each block in a specific 
order, starting with the secondary blocks at the very end of the decay chain, and ending 
with the main block. 

This procedure is best illustrated with the example in Figure |6|. A phase-space point 
is defined by generating all the integration variables in the transformed expression of the 
phase-space measure: the invariant mass of each leg shown as a dashed line, the direction 
{0, (p) of any visible particle, and the energy of the visible particles represented by the solid 
thin lines. Then all other kinematic variables are determined as a function of the generated 
variables, first in the secondary blocks A and E, then in the secondary block D (by means 
of formulas that are parametrized by the kinematics in block E), and finally in the main 
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block B (by means of formulas that are parametrized by the kinematics of all the secondary 
blocks). Such an approach can be easily generalized to the case of an arbitrary decay chain. 
One difficulty in our approach is that the boundary of the physical phase-space vol- 
ume cannot be translated into simple conditions on the variables of integration. In order 
to implement the boundary, we simply check point-by-point in the phase-space that the 
variables in the main block are physical (for example, the Bjorken fractions qi or q2 can- 
not be larger than one), otherwise we throw away the phase-space point. In some cases, 
the fraction of unphysical points that are removed in this way may be large. Still, the 
algorithm is rather fast since the generation of a phase-space point is in general much less 
time-consuming than the evaluation of the squared matrix element. 

3.4 Mult i- channel phase-space generator 

If a given parametrization of the phase-space measure maps all the peaks in the integrand 
simultaneously, an adaptive Monte Carlo integration using only this channel is expected 
to be efficient. But most of the time, each peak in the integrand cannot be mapped onto a 
variable of integration in a single phase-space mapping, since the number of peaks is larger 
than the dimension of the phase space. ^ In these cases, we keep several channels, i.e. 
several phase-space parametrizations z ^ z' = Pi{z) such that each peak in the integrand 
is mapped onto a variable of integration in at least one channel. The total integration 
can be carried out using a multi-channel integration approach in which every channel i 
comes with a phase-space-dependent weight Pi{z) > in the global parametrization of the 
phase-space measure: 

z^z' = P{z) = Y,Pi{z)Pi{z), (3.5) 

i 

with the condition ^^^Piiz) = 1. Each weight /3i(z) must be chosen such that it is signifi- 
cant in the phase-space region where the corresponding channel Pi{z) is relevant. In the 
case of the computation of total cross sections, this condition can be automatically fulfilled 
by setting I3i{z) to be proportional to the amplitude squared of a single diagram associated 
with the channel i [?]. In analogy to the single-diagram enhanced method, we choose to 
set the weight I3i{z) to be proportional to the product of the peaks that are mapped onto 
integration variables in the corresponding phase-space mapping. 

In comparison with previous implementations for the evaluation of the matrix element 
weights, this multi-channel approach is expected to speed up the convergence of the in- 
tegration, especially in the case of an over-constrained topology. An example of such a 
topology has been investigated in [?, ?, ?], where either the helicity of the W boson or the 
mass of the top quark is reconstructed from ti events in the semi-leptonic channel. In these 
analyses, a single channel was used for the evaluation of the weights, leaving unmapped 
a subset of peaks in the integrand. On the contrary, our procedure always maps a given 
peak in the integrand onto a variable of integration in at least one channel. 



^This situation could also occur if one of the required blocks to build such a phase-space mapping 
corresponds to a change of variables that cannot be inverted analytically and hence has not been considered 
in our algorithm. 
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blocks 


integrated volume 


3 


MB A 


6.30 X 10-5 


3 


MB B 


6.30 X 10-5 


3 


MB C 


6.30 X 10-5 


6 


MB D 


694 GeV^ 


4 


MBE 


0.0166 GeV^ 


4 


MBF 


0.0166 GeV2 


5 


MB B + SB A 


3.89 GeV^ 


4 


MB B + SB B 


0.0166 GeV2 


3 


MB B + SB C 


6.30 X 10-5 


3 


MB B + SB D 


6.30 X 10-5 


4 


MB B + SB E 


0.0166 GeV^ 



Table 1: Phase-space volumes ^ dqidq2d(j)nl / {sqiq2) for I massless particles produced in hadron- 
hadron collisions at y^ = 1 TeV. The number I of final-state particles is indicated in the first 
column. The second and third columns indicate the structure in blocks defining the phase-space 
mapping that is used to calculate the volume with our phase-space generator, and the numerical 
value that wc obtained at three digits accuracy. We check that those number agree exactly with 
the exact value of the phase-space volume. 

The whole procedure that we have presented so far has been implemented in the 
MadGraph framework, and the corresponding module has been named MadWeight. For a 
given decay chain and a transfer function for the final state objects, the optimized phase- 
space mappings are automatically selected, and the resulting multi-channel phase-space 
generator is used for the evaluation of the weights. While this procedure applies for virtually 
all cases, the speed of convergence of the numerical integration strongly depends on the 
process under investigation, and whether the calculation time is a serious limitation or not 
has to be assessed on a case- by-case basis. 

3.5 Validation of the phase-space generator 

One potential issue related to our phase-space mappings optimized for the computation of 
the weights is the fact that some of the associated jacobians develop singularities in specific 
phase-space regions. These singular regions are an artefact of the change of variables. In 
our case they have a null measure in the integration volume. One can therefore split the 
integration volume into a volume Vi where the jacobian is finite and a volume V2 that 
contains the singular region and that can be made arbitrary small compared to the volume 
Vi. At any given accuracy, we can ignore the contribution from the volume V2 provided 
that e = V2/V1 is sufficiently small. At the numerical level though, one may fear that 
instabilities will appear in this procedure. 

In practice, we have not encountered any numerical instabilities resulting from a change 
of variables that is associated with a specific phase-space block. Any phase-space block and 
the related change of variables that have been defined in our procedure have been checked 
by reproducing the volume of the entire phase-space region with our phase-space generator 
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I 


blocks 


integrated volume 


3 


MB A 


3.49 X 10-5 


3 


MB B 


3.49 X 10-5 


3* 


MB C 


4.13 X 10-5 


6 


MB D 


124 GeV^ 


4 


MBE 


8.17 X 10-3 GeV^ 


4 


MBF 


8.17 X 10-3 GeV^ 


5 


MB B + SB A 


1.28 GeV^ 


4 


MB B + SB B 


8.17 X 10-3 GeV^ 


3 


MB B + SB C 


3.49 X 10-5 


3 


MB B + SB D 


3.49 X 10-5 


4* 


MB B + SB E 


9.78 GeV^ 



Table 2: Phase-space volumes J dqidq2d(j)nl / {sqiq2) for I particles with a mass ?ti = 50 GeV 
produced in hadron-hadron collisions at y/s — 1 TeV. The number / of final-state particles is 
indicated in the first column. A star * indicates that the mass of one of the final state particles is 
set to zero, as this condition is required by one of the blocks. The second and third columns indicate 
the structure in blocks defining the phase-space mapping that is used to calculate the volume with 
our phase-space generator, and the numerical value that we obtained at three digits accuracy. We 
check that those number agree exactly with the exact value of the phase-space volume. 

using a parametrization of the phase-space measure that involves this block. This Monte 
Carlo procedure to compute the phase-space volume has a very poor convergence, as the 
phase-space mappings that are optimized for the computation of the weights are clearly 
inefficient for the computation of just the phase-space volume. Nevertheless, by increasing 
the number of generated phase-space points, we checked that the phase-space volume is 
reproduced with an accuracy better than one percent for each tested phase-space mapping. 
We first set the mass of the final-state particles to zero and obtained the results summarized 
in Table ||. We then considered the case of massive particles in the final state and obtained 
the results summarized in Table |2[ 

In order to validate the multichannel implementation, we also computed the total cross 
section of several processes by integrating the squared matrix element with our phase-space 
generator. This can be achieved by setting all transfer functions to one. Here again, the 
convergence of the numerical integration is poor, as the phase-space parametrization is not 
designed for such computation. By using a very high statistics, we reproduced the total 
cross sections associated with the processes listed in the first column of Table |3|. 

4. Example of applications 

In this section we illustrate a few examples of studies that can be achieved with Mad- 
Weight. The following analyses are based on simulated events generated with Mad- 
Graph/MadEvent [?]. The events are passed through Pythia [?] for the showering and 
the hadronization. Electrons and muons are assumed to be reconstructed with 100% ef- 
ficiency and with excellent resolution if they have a pseudo-rapidity \r]\ < 2.4. Detector 
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process 


^MW/^ME 


channels 


blocks 


pp-^{W ^ jj)j 


0.982(6) 


3 


MBA 


pp-^ {W ^ Iv) 


0.9991(14) 


1 


MBB 


pp^ [W ^9r{f- >T-x)] 


1.003(5) 


1 


MB B; SB C 


pp -^ 2[fi -^ iix] 


1.020(5) 


3 


MB B,F; SB C 


pp -^ 2[t -^ b{W -^ li^i)] 


1.000(25) 


1 


MBD 


pp^[t^ b{W+ -^ lui)] + [i^ b{W- -^ jj)] 


0.94(5) 


6 


MB B, SB D,E 


pp^h^ {W+ -^ fl+Um){W- -^ 11- Urn) 


0.99(2) 


1 


MBE 



Table 3: Validation of the phase-space generator by computing total cross sections. The processes 
under consideration are written in the first column. The second column gives the ratio of the cross 
section computed with MadWeight over the one computed with MadEvent [?]. The third column 
indicates the number of channels that are used in the MadWeight integration, and the last column 
indicates the blocks that are involved in that integration. 



response simulation is performed using PGS [?] which takes into account geometrical ac- 
ceptance, finite granularity and energy resolutions of typical calorimeters used in LHC 
experiments. Jets are then reconstructed based on the kt algorithm [?, ?, ?] and applied 
on the calorimeter cells fired by the generated stable or quasi-stable particles. 

The transfer functions VFj (x*,y*), Eq. ( p.4D with i running over all reconstructed jets 
are determined from an independent ii sample where well separated jets (including light and 
b jets) are matched to the corresponding partons. We consider a double- Gaussian shape 
function characterized by 5 parameters: the means and the widths of the two Gaussian 
distributions, and their relative normalization. We fit these five parameters in each 20 
GeV bin in jet energy from 40 GeV to 200 GeV. The energy dependence of the mean 
and the width of each Gaussian distribution is then approximated by the parametrization 
Cl+C2^fE+c^E, with the coefficients ci, C2 and C3 extracted from a x^ ^t to the values of the 
four parameters of Gaussian distributions in each energy bin. The relative normalization 
of the two Gaussian distributions is assumed to be energy independent, and is fixed to 
the average of the corresponding values in each energy bin. The typical resolution for jet 
energy is between 5 and 12 GeV, with tails parametrized by Gaussian of variances as large 
as 30 GeV. 



4.1 Top-quark mass measurement 

The top-quark mass measurement by means of the matrix element method was published 
for the first time by the D0 collaboration using the single-leptonic final state arising from 
top-quark pair production [?]. The method has been later extended also to include a 
simultaneous determination of the Jet-Energy-Scale uncertainty. The accuracy of the ex- 
perimental determination of mt has been further improved by the contribution of other 
studies based on matrix element method [?,?,?]. In these analyses, a dedicated phase- 
space integration was performed to define the event weight. Our algorithm provides this 
weight automatically based on the blocks given in Table |3|, last but one line. 

As an example of application of our automatic reweighting algorithm, we illustrate the 
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Figure 7: (a) Logarithmic likelihood values for a sample of 20 events generated with rrit^ input = 170 
GeV. The solid line is a parabolic fit to the points near the minimum. The statistic error is 
estimated by the half width of the distribution at log(L/imax) = 0.5 and is extracted from the 
fit. (b) Calibration of the matrix element mass fitting procedure. The errorbars correspond to the 
the value of mass of the top quark and the associated statistic error reconstructed from tt samples 
generated with different input values of rrit- The solid line is a linear fit to the four points and the 
dotted line corresponds to mice = f^input- 

performance of the method for the determination of the top-quark mass at the LHC, by 
using a small statistics of tt events in the single lepton final state: 



p'p 



[i^ b{W- ^ t,-u^][t ^ biW+ ^ jj)]. 



(4.1) 



For the sake of simplicity, we assume that there is no background and 20 signal events 
after selection. Pseudo-data have been simulated with an input top-quark mass at 170 
GeV. The selection requires one muon with a reconstructed transverse momentum above 
10 GeV and exactly four isolated jets with a reconstructed transverse momentum above 20 
GeV. 

The determination of the top-quark mass from our sample of pseudo-data is obtained 
by the minimization of — log(L) with respect of rrit where the likelihood L is defined -up 
to a normalization factor- by the product of the weights calculated for each event 



N 



log(^) = -^ log[P{xi;mt)] . 



(4.2) 



i=l 



The acceptance of the detector and the cuts imposed on the sample may depend on the 
input mass of the top quark. Such a dependence might introduce a bias in the extraction 
of m-top from the fit of the likelihood given in Eq. ([4.2|). We explicitly tested that in our 
pseudo-data this bias is very small and we therefore ignored it in this example. 

The values of — log[L(?ni)] for different assumptions of mt are displayed in Figure |^( a). 
A clear minimum is observed close to the input mass value. A parabolic fit gives the 
value rrit = 171.9 ± 2.0stat GeV. Using ten independent samples generated under the same 
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conditions, we extracted the expected mt value and the expected statistic error for such a 
measurement: mt = 173.5 ±3.7 stat GeV. We also identified a small significant shift between 
the input and the reconstructed mass by repeating the analysis for other input masses, as 
can be seen from the calibration curve displayed in Figure |^(b) . This curve was obtained by 
generating tt samples with top quark masses of 160, 170, 180 and 190 GeV with the same 
selection procedure and the same fitting procedure, and with a statistics of 100 events per 
sample. The small bias between the input and the reconstructed mass of the top quark can 
possibly result from the effect of initial state radiation, which is not completely removed 
in our event selection criteria [?]. 

4.2 Spin identification in decay chains with missing energy 

As a second illustration we address the challenge of determining the spin of new particles. 
A simple example is the production of a light charged Higgs boson with a mass close to the 
W boson mass. The information on the spin of a resonance is passed through the angular 
distribution of its decay products. If the momentum of each final-state particle produced 
in the decay chain is measured, the angular distributions can be reconstructed and the 
spin of the resonance identified. In fact in many cases, such as the one we have chosen, 
the final state is characterized by missing transverse energy from undetected particles 
and the angular distributions of the decay products cannot be fully determined. The 
interesting question becomes therefore whether the available information from the final 
state is sufficient to discriminate between different spin assignments in the decay chain. 
The matrix element method appears to be particularly relevant in this case, since the event 
weight will encompass the whole available event kinematics including the spin correlation 
effects that survive after the experimental reconstruction of the events. 

In the following example, we assume that the production of the signal and its irreducible 
background proceeds exclusively via the production of top quark pair that subsequently 
decay into H~^ + 6 or into W + b, with mfj± ~ m^r± , followed by a leptonic decay of both 
bosons. The signal process is 

pp^[t^ b{H+ -^ T+Vr)\[i^ b{W- -^ ^l-u^,)l (4.3) 

and the corresponding irreducible background results from the production of a pair of W 
bosons 

pp^[t^ b{W^ -^ T+Ur)][i^ b{W- -^ ^i'i>^,)]. (4.4) 

At the reconstruction level, we required the presence of exactly two jets with a pT larger 
than 20 GeV, one r"*" — assumed to be reconstructed as precisely as the other charged 
leptons — and one /x~. The cuts on these leptons are |ry| < 2.4 and px > 5 GeV. We 
reject the events containing photons or electrons with pj- > 5 GeV and |?7| < 2.4, while 
no restriction is imposed on the number of jets with a px less then 20 GeV. The transfer 
functions associated with the other particles have the same parametrization as described in 



Section 4.1. We arbitrarily choose a final relative normalization of signal and background 



events, working with a sample of 240 signal events and 760 background events. 
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Figure 8: Expected normalized distribution of events with respect to the discriminant d built upon 
(a) the matrix element weight and (b) the px of the tau for a pure signal sample (solid histogram) 
and for a pure background sample (dashed histogram). The errorbars are the distribution associated 
with the pseudo-experiment sample, assuming that the statistical error on the number N of events 
in a bin is given by VN. 



Before proceeding further, we stress that while providing an interesting case study, our 
example cannot be regarded very realistic. First, the relative cross sections and reconstruc- 
tion efficiencies for the signal and the background have been chosen arbitrarily. Second, 
this such a light charged Higgs is not favoured by the present constraints, which point to 
much higher masses. Finally, the tau lepton reconstruction is idealized as it is considered 
here on the same footing as as a muon. A more realistic approach would consist of taking 
into account the energy loss from the tau decay with a dedicated transfer function for the 
energy of the tau. Nonetheless, as shown below, this example illustrates quite well the 
power of the matrix element method. 

Let us define Psix), Pb{x) as the weights evaluated for the event final state x under 
the signal and the background hypotheses, respectively. These weights can be calculated 
from the signal and background full matrix element as defined in Eq. ( |2.2| ). Alternatively, 
they can be associated with a normalized differential cross section with respect to a single 
observable, such as the r"*" transverse momentum 



Ps,b{x) 



1 das,B 



[PTir"-)] 



(4.5) 



(^S,B dpT 

which also captures the spin effects. The advantage of the weights defined in Eq. (|4.5|) is 



in their simplicity: their evaluation only requires to use a standard phase-space generator 
that is optimized for the computation of cross sections. Such an observable, for example, 
is very commonly used in the determination of the polarization of the W bosons in top 
events and provides us with a useful benchmark to study the increased sensitivity that the 
matrix element method might provide. 

These weights can then be combined to build an event-by-event discriminating variable 

Ps{x) 



d{x) 



Ps{x) + Pb{x) 



(4.6) 
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Figure 9: x^ values associated with the fit of the pseudo-experiment data to the theoretical 
prediction parametrized by the fraction r of signal events for the weights calculated from (a) the 
matrix element, (b) the pt of the tau. Each dashed line represents the one-standard-deviation 
interval defined by the condition x^(r) < x^j^^ -I- 1. 



The normalized-to-one distributions of events as a function of the discriminant variable d 
are shown in Figure ^ for the two cases. The solid (resp. dashed) histogram is the distri- 
bution expected for a pure sample of signal (resp. background) events. These distributions 
have been generated from large samples, in order to allow us to neglect the statistical 
fluctuations. Spin correlation effects are expected to give rise to different values for the 
weights under the two spin hypotheses. Nevertheless, for most of the events, this disparity 
is expected to be small, resulting in a discriminant close to d ~ 0.5. Note that in our 
example, signal and background events are characterized by the same topology with inter- 
mediate particles of the same mass. Only the spin of the intermediate W or H resonances 
differ between the two decay chains. The distributions clearly shows that the discriminant 
power is substantially reduced when only the information on the transverse momentum of 
the T~^ is retained. 

Yet, for some events the discriminant is significantly different of 0.5, corresponding to 
configurations clearly favoured by one of the two hypotheses. Such events influence the 
shape of the distributions and allow us to distinguish them. One can take advantage of this 
difference to find out the fraction of signal events in the pseudo-experiment sample. The 
normalized-to-one distribution associated with the pseudo-experiment sample as a function 
of the discriminant variable d is also displayed in Figure ^. The fraction of signal events in 
the pseudo-experiment sample can be reconstructed by a least-square fit, i.e. by minimizing 



x'(r) 



bins in d 



{-Pdatajd) - [rVsjd) + (1 - r)VBid)]y 

[AVdataidW 



(4.7) 



where Vs, Vb, are the expected binned distributions for signal and background events and 
Vdata is the binned distribution associated with the pseudo-experiment sample. 

The x^ values as a function of the fraction of signal events are shown in Figure p. The 
best fit for each discriminant are obtained for r = 24 it 9% and r = 30 it 23%, respectively. 
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Figure 10: Generic decay chain corresponding to the production of two resonances followed by 
their decay into weakly-interacting and standard model particles. The weakly-interacting particles 
are represented by the dashed lines. 

Both results are compatible with the true fraction of signal events, but the discriminant 
using matrix elements increases the accuracy by more than a factor of two. 

4.3 Smuon pair production at the LHC 

Over the past fifteen years, a tremendous amount of work has been devoted to new tech- 
niques for mass reconstruction of new particles that might be produced at the LHC. Accord- 
ing to most scenarios, the hypothetical new physics states are not expected to be directly 
observed experimentally, i.e., they appear as intermediate states in specific decay chains or 
they escape from the detector without interacting with it. Their mass can hence only be 
reconstructed indirectly, by making a number of assumptions on the decay chain at work. 
The number of assumptions in turn is directly correlated to the amount of information 
that can be extracted from the decay chain. However, due to the lack of constraints on 
physics beyond the standard model, the proposed techniques have to be general enough, 
at least if they are aimed at reconstructing the mass of new hypothetical particles in the 
early stages of investigation. Further more, the limited knowledge of the detector has to 
be taken into account. It is in this context that a number of mass measurement techniques 
based on kinematic methods have been proposed in the literature. They can be classified 
according to the type of decay chains that they address and according to the assumptions 
on which they rely [?]. 

Despite the plurality of kinematic variables that have been proposed, mass determi- 
nation remains very challenging for specific decay chains. One well-known example of a 
difficult topology is the production of two resonances followed by their decay into a weakly- 
interacting and a standard model particles, shown in Figure |l^. In this case, the kinematic 
methods that have been proposed to reconstruct simultaneously the mass of the two new 
particles require a very high statistics. Whether their sensitivity is sufficient under real 
experimental conditions still remains to be determined. A complementary way to address 
the same problem is to ask what would be the maximum sensitivity achievable, given a 
very detailed set of hypothesis to be tested that not only include masses but also the spin 
and couplings information, i.e., taking into account the full theoretical model prediction. 

The problem can be investigated with the matrix element method that usually makes 
use of the strongest assumptions on the analysed events [?]. One way to dramatically 
increase the theoretical information is to assume that the masses of the new physics states 
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Figure 11: Logarithmic likelihood as a function of the hypothesis values for [(tti- — 
TO? )/2'mjx^,my--^] built upon (a) the matrix element weights, (b) the transverse momentum of 
the fi^ and the invariant mass of the muons. 

are the only unknown properties of the decay chain. We therefore consider a specific decay 
chain corresponding to the topology in Figure |l^: the production of a pair of smuons 
followed by their decay into a muon and a neutralino 

pp^ (/i+ ^^+Xi)(Ar ^M^Xi)- (4.8) 

We suppose that we have isolated a pure sample of events that correspond to the decay 



chain in Eq. (4.8). Further, we assume a perfect reconstruction of the kinematics of the 
two muons in each event. Within these assumptions, the significance that can be achieved 
with the matrix element method provides us with an upper bound on the significance that 
can be delivered by any realistic analyses at a given luminosity. 

We have considered the following input hypothesis for the masses of the sparticles: 

m-/i,,mput = 150 GeV, m^,, input = 100 GeV. (4.9) 

Under this hypothesis, we have generated events corresponding to the decay chain in 
Eq. ([4.^). We have built a sample of fifty events with exactly one fi'^ and one ^~ with a 
transverse momentum larger than 20 GeV, and no other particles except maybe some jets 
with a pt less than 20 GeV. These events are regarded as a pseudo-experiment sample in 
the following. 

The sensitivity that can be achieved with the matrix element method has been anal- 
ysed by computing the weights P(xj|m^^,?TT,^^) for each event Xi in the pseudo-experiment 
sample. A bias may be introduced by the acceptance cuts. This effect has been corrected 
by normalizing the probability density in the acceptance region. This amounts to replace 
the factor 1/aa by the factor l/a^^ in the definition of the probability density in Eq (^), 
with a the cross section in the acceptance region. In terms of the weights normalized in 
the acceptance region, the unbiased likelihood has the usual form 

iV=50 



log L{mf,^ 



y^ logP(xj|mp 



(4.10) 
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It is advantageous to express the likelihood in term of the variable (m? — m3- )/2m^^, 
that corresponds to the momentum of each final state particle in the rest frame of the 
smuon from which it originates. The complementary variable can be chosen to be "i^^. 
The likelihood for different theoretical hypotheses is shown in Figure |ll|(a). The optimal 
value for the variable (fn? — m|^)/2m,^^ is 42 GeV, which corresponds to the input value. 
There is a very mild sensitivity with respect to variation of the complementary variable. 

One way to highlight the increase of sensitivity by using the complete theoretical and 
experimental information is to compare the profile of the likelihood built upon the matrix 
element weights (shown in Figure 0(a)) with the likelihood profile that is obtained by 
keeping only the information contained in the transverse momentum px^ of the /i"*" and 
the invariant mass M^^ of the muons. In order to simplify the computation of this second 
likelihood profile we neglect the correlations between the two variables pxfi and M^^. Thus 
the weight attached to each event is reduced to 

, I , 1 da . . \ 1 '^'^ / 71 J- I \ /A -,-,\ 

The resulting likelihood profile is displayed in Figure pT|(b). The comparison with Fig- 
ure |ll|(a) shows that the sensitivity to the theoretical hypothesis (m^^ , rrij^^ ) is dramatically 
reduced when only the information contained in the kinematic variables pT^ and M^^ is 
used. 

5. Conclusion 

We have presented a new algorithm that allows the automatic computation of the weights 
appearing in the matrix element method. Given an arbitrary decay chain and a transfer 
function tuned to the resolution of the detector, our code produces a specific phase-space 
generator that combines different phase-space mappings optimized for the integration of 
the product of the matrix element and the transfer function. The mappings are obtained 
by applying a specific transformation to the standard parametrization of the phase-space 
measure. This transformation is expressed in terms of a composition of changes of variables 
acting on different kinematic sectors of the topology. As a result, our algorithm leads to 
a modular structure, and hence is particularly convenient for future improvements. For 
example, the current implementation could be easily extended to include non-analytical 
changes of variables. 

The availability of a tool such as MadWeight that provides the resource for the auto- 
matic evaluation of the weights, sparing the user to focus on the technical details of matrix 
element generation and integration over phase space, paves the way to a potentially large 
number of new applications. First MadWeight could be used to improve our understanding 
of the matrix element method itself and of its limitations. For example, its implementation 
within Madgraph provides all the required computational tools to analyse the infiuence 
of additional jet radiation in a specific measurement, or to estimate the systematic error 
resulting from the parametrization of the transfer function. Second not only measurement 
of masses or cross sections in Standard Model could be achieved in a effortless and more 
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efficient way, but the matrix element method could also be employed in the search and 
identification of new physics models. Different hypotheses, such as those corresponding to 
(any) new physics scenario and/or benchmark parameter points could be tested against 
data and be assigned a meaningful relative probability. We look forward to exciting new 
developments in these directions. 
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A. Phase-space measure associated with the main blocks 

A.l MB A 

The notation for the phase-space variables associated with this main block is given in 



Figure 12. The two momenta pi and p2 correspond to the visible particles that enter into 

'?! -^ / / bil 
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Figure 12: Notation for the kinematics of MB A. 

the main block, along with the Bjorken fractions qi and 52- The standard phase-space 
parametrization associated with this main block reads 

The four-vectors Pjn and Pgn refer to the total momenta in the initial and final states, 
respectively. In our procedure we apply a change of variables that leads to the following 
parametrization of the phase-space measure 

167r2^ E dMhde2d^2 x J, (A.2) 



27 



where 9i and (j)i refer to the polar and azimuthal angles of particle i with respect to the 
beam axis. The jacobian J of this transformation reads 



T "^ I i2i i2i 

J = -|Pir|P2rl cos( 

s 



)i sm (p2 — sm (pi COS (p2 1 



(A.3) 



where s is the squared invariant mass of the colliding hadrons. The energies Ei, E2 of 
the final particles in the main block are adjusted to balance the transverse momentum 



pbranchcs q£ ^t^ ^]^g branches represented by the blobs in Figure 12. We assume that this 



transverse momentum is different from zero, except maybe in a region of null measure (see 
the discussion in Section |3^ ) . This requires the number of particles in the final states to be 
larger than or equal to three. Then the variables \pi\, IP2I ^^"^ be expressed as the solution 
of the following linear system 



\Pi I sin 61 cos (pi + |p2 1 sin O2 cos <p2 
\Pi I sin 9i sin (f)i + |p2 1 sin 62 sin (f)2 



_ branches 
Px 1 



-p, 



branches 

y 



(A.4a) 
(A.4b) 



The Bjorken fractions qi,q2 are then fixed by imposing the conservation of total energy 
and total momentum along the beam axis. 

A.2 MB B 

The notation for the phase-space variables associated with this main block is given in 



Figure 13. The momentum pi corresponds to the missing particle that belongs to the main 
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Figure 13: Notation for the kinematics of MB B. 

block, along with the Bjorken fractions qi and (72- The momentum p2 corresponds to the 
branch that is directly connected to that missing particle. The variable S12 is the invariant 
(pi + ^2)^- The standard phase-space parametrization associated with this constrained 
sector reads 



dqidq2 



-(27r)V(Pi„-Pfin) 



(A.5) 



' (27r)32^i 

The four-vectors Pjn and Pgn refer to the total momenta in the initial and final states, 
respectively. In our procedure we apply a change of variables that leads to the following 
parametrization of the phase-space measure 



A-kEx 



dsi2 X J. 



(A.6) 
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The jacobian J of this transformation is given by 



J = 1^22-^1 - E2PIZ 



(A.7) 



where s is the squared invariant mass of the colhding hadrons. The transverse momentum 
the missing particle is fixed by requiring that it balances the transverse momentum of all 



the branches represented by the blobs in Figure |13|. The component piz of momentum 
along the beam axis is fixed by imposing the invariant mass condition 

{Pl+P2f 



Sl2- 



(A. 



If the energy Ei of the missing particle is treated as an independent parameter, the left side 
of the Eq. ( |A.8| ) is a first-order polynomial in pi^. We therefore obtain a unique expression 
for piz in term of Ei. The mass-shell condition associated with the missing particle gives 
rise to up to two solutions for the energy E\. Each solution that gives a real positive value 
for El and that leads to values of the Bjorken fractions qi and §2 between and 1 is kept, as 



it corresponds to a distinct physical phase-space point at which the Jacobian in Eq. ([A.?]) 
and the integrand must be evaluated. 

A.3 MB C 

The notation for the phase-space variables associated with this main blob is given in Fig- 



ure 



14. The momentum of the missing particle is denoted by pi, the momentum of the 
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Figure 14: Notation for the kinematics of MB C. 

branch directly connected to the missing particle is denoted by p2^ the momentum of the 
massless visible particle in the main block is denoted by ^3. The Bjorken fractions are 
denoted by qi and 52- The variables S12 and S123 refer to the invariants {pi + ^2)^ and 
[pi + P2 +^3)^1 respectively. The standard phase-space parametrization associated with 
this main block reads 

d^pi d?Pz 



dqidq2 



(27r)^5* (Pin - Pfin) 



(A.9) 



' (27r)32£;i i2TT)^2E3 

The four-vectors Pjn and Pgn refer to the total momenta in the initial and final states, 
respectively. In our procedure we apply a change of variables that leads to the following 
parametrization of the phase-space measure 



16tt'^EiE3 



id9sdsi2dsi23 X J- 



(A.IO) 
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The jacobian J of this transformation is given by 

XE2PIZ - xEiV2 



J = sin03 



^l^i 



s 
2cos{4>^) cos 



, ^ , (03)^2^1x^3 sin(03) + 2 cos(03) cos(03)^iP2x£^3 sin(e3) - 

2cos{<i)z)pizP2xE-i sin(6'3) + 2 (zos{4>:>,)pixP2zE:i s\n{6^) - 

2 cos{e3)E2PiyE3 sin((/)3) sin(6'3) + 2 cos{9:>,)Eip2yE-i sin((/)3) sin(6'3) - 

2pizP2yE3 sm{(j)3) sin(6'3) + 2piyP2zE3 sin((/)3) sin(6'3) + 

2cos((/.3)2^2Pi.^3sin(e3)^ - 2cos(</.3)2£;ip2^^3sin(03)^ + 

-1 

2E2P1ZE3 sin(</.3)2 sin(03)^ - 2Eip2zE3 sin(</)3)2 sin(03; 



(A.ll) 



with X = 2p3.(pi + P2)/E3 and s standing for the squared invariant mass of the colhding 
hadrons. 

If we treat the variables Ei and a = 2pi.p3 as two independent parameters, the 
components of the three momentmn Pi of the missing particle and the energy E3 = \p^\ of 
the massless visible particle can be expressed as the solution of the following linear system 
of four equations 

{Pi+P2? = si2 (A.12a) 

{Pl + P2 + PS? = Sl23 (A.12b) 

pi, + E3 sin 03 cos </.3 = -p'^r'"'' (A. 12c) 

Ply + E3 sin 03 sin ^3 = -P^'^''''''' (A. 12d) 

that is parametrized by the momentum p2, by the angles ^3 and (p3, by the total transverse 
momentum ^branches q£ g^ij ^]^g branches represented, by the blobs in Figure 14 and by 



the variables a and Ei. The next step is to determine the values of the variables a and 
El. The mass-shell condition for the missing particle of momentum pi and the equation 
2pi-P3 = ct defines a system of two coupled quadratic equations in the variables Ei and 
a, parametrized by the momenta of the blocks. This system can be solved analytically. 
There are up to four solutions for Ei and a. Each solution that is physical {i.e., such that 
\p3\ > 0, El > and each of the Bjorken fractions qi,q2 is between and 1) corresponds to 
a distinct phase-space point at which the Jacobian in Eq. ( |A.ll ) and the integrand must 
be evaluated. 

A.4 MB D 

The notation for the phase-space variables associated with this constrained sector is given 



in Figure 15. The momenta of the missing particles are denoted by pi and p2, the momenta 
of the branches connected to the main block are denoted by p3, p4,, p^ and pQ. The Bjorken 
fractions are denoted by qi and 52- The variables Sij and Sjjfc refer to the invariants {pi+Pj)"^ 
and (pi+Pj+Pk)'^, respectively. The standard phase-space parametrization associated with 
this main block reads 
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Figure 15: Notation for the kinematics of MB D. 



The four-vectors Pjn and Pgn refer to the total momenta in the initial and final states, 
respectively. In our procedure we apply a change of variables that leads to the following 
parametrization of the phase-space measure 



The jacobian J of this transformation is given by 

E1E2 



dsi3dSi34,dS25dS256 X J- 



J 



8s 



E3IE5 [p3Az{PlyP2zP5ex " PlxP2zP56y 



-PlyP2xPmz +PlxP2yP56z) + Plz{-P2zP34yP56x + 
P2zP34xP56y - P2yP34xP56z + P2xP3AyP56z)] + 
{E^6P2z - E2P56z){PlzP34yP5x " PlyP34zP5x " PlzP34xP5y + 
PlxP3AzP5y) + [E56{PlzP2yP34x - PlzP2xP3Ay + PlyP2xP34z " 
PlxP2yP34z) + E2{pizP3AyP56x " PlyP34zP56x " PlzP34xP56y + 
PlxP3AzP56y)]P5zj + EsAE5P2z(pizP3yP56x - PlyP3zP5Gx 
-PlzP3xP56y +PlxP3zP56y) + E^{pizP2yP3x - PlzP2xP3y 
+PlyP2xP3z - PlxP2yP3z)P56z - {E^&P2z - E2P56z) 
iPlzP3yP5x - PlyP3zP5x - PlzP3xP5y + PlxP3zP5y) 
-[E56{PlzP2yP3x - PlzP2xP3y + PlyP2xP3z - PlxP2yP3z) + 
E2ipizP3yP56x - PlyP3zP56x - PlzP3xP56y + PlxP3zP56y)]p5z> 

El\[E5{p2z{-p3AzP3yP56x + P3AyP3zP56x + 
P3AzP3xP56y - P3AxP3zP56y) + 

{-P2yP3AzP3x + P2xP3AzP3y + P2yP3AxP3z - P2xP3AyP3z)P56z] + 
[E56P2Z - E2P5(iz){P3AzP3yP5x - P3AyP3zP5x - PSAzPSxPby + 
P3AxP3zP5y) + {E^6{P2yP3AzP3x - P2xP3AzP3y - P2yP3AxP3z + 
P2xP3AyP3z) + E2{p3AzP3yPmx - P3AyP3zP56x - P3AzP3xP56y + 



+ 



P3AxP3zP56y 



)]P5z} 



(A.14) 



(A.15) 
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where Eij = Ei + Ej , pij = pi + pj , and s is the squared invariant mass of the colhding 
hadrons in their center-of- momentum frame. If we treat the variables Ei and E2 as in- 
dependent parameters, then the components of the three momenta Pi,P2 of the missing 
particles can be expressed as the solution of the following linear system of six equations 



{pi + psY' 



Sl3 



(Pl +P3+P4.f = S134 

{P2+P5f = S25 
{P2+P5+P6f 



•S256 



^ _L ,r, — ^branches 
Plx + P2x — -PTx 



Ply + P2y 



^branches 
-PTy 



(A.16a) 
(A.16b) 
(A.16c) 
(A.16d) 
(A.16e) 
(A.16f) 



that is parametrized by the momenta p^, . . . ,pQ of the branches connected to the main 
block, by the total transverse momentum p™-anches q£ ^y[ the branches represented by the 



blobs in Figure 15 and by the variables Ei and £'2- The next step is to determine the 
values of the variables Ei and E2 ■ The mass-shell conditions for the two missing particles 
of momentum pi and p2 define a system of two coupled quadratic equations in the variables 
El and E2, that can be solved analytically. There are up to four solutions for Ei and £^2- 
Each solution that is physical (i.e., such that E2 > 0, Ei > and each of the Bjorken 
fractions qi,q2 is between and 1) corresponds to a distinct phase-space point at which 
the Jacobian in Eq. ( |A.15 ) and the integrand must be evaluated. 



A.5 MB E 

The notation for the phase-space variables associated with this main blob is given in Fig- 



ure 16. The momenta of the missing particles are denoted by pi and p2, the momenta 





Figure 16: Notation for the kinematics of MB E. 

of the branches directly connected to these missing particles are denoted by p^ and p4. 
The Bjorken fractions are denoted by gi and q2- The variables Sij refer to the invariants 
{pi +Pj)'^ and s denotes the squared invariant mass of the colliding partons. The standard 
phase-space parametrization associated with this MB reads 



dqidq2 



cP'pi (fp2 

(27r)32£i (27r)32£2 



(27r)^5^ (Pin - Pi 



fin J 



(A.17) 
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The four-vectors Pin and Pgn refer to the total momenta in the initial and final states, 
respectively. In our procedure we apply a change of variables that leads to the following 
parametrization of the phase-space measure 

-dydsdsi^ds2i x J, (A. 18) 



where y is the rapidity of the colliding partons in the lab frame. The jacobian J of this 
transformation is given by 

E1E2 
" 4s 

PlxP2yP3z) + E2PlzP3yPix " EiP2zP3yPix " E2PlyP3zP4.x + EiP2yP3zPix 

-E2PlzP3xP4y + EiP2zP3xP4y + E2PlxP3zPAy " EiP2xP3zP4y + (-E'2Plj/P3x + 

-EiP2yP3x - E2PlxP3y + EiP2xP3y)Piz + E2,{-pizP2yPix + PlyPlzPix 



J 



EA{PlzP2yP3x - PlyP2zP3x " PlzP2xP3y + PlxP2zP3y + PlyP2xP3z 



+PlzP2xP4y - PlxP2zPAy - PlyP2xPAz + PlxP2yPAz 



(A.19) 



If we treat the variables £'1, E2 and p2y as independent parameters, the other compo- 
nents of the momenta pi , p2 of the missing particles can be expressed as the solution of the 
following linear system of five equations 

{pi+P3? = si^ (A.20a) 

{P2+PA? = S2A (A.20b) 

Pix+P2x = -p''J''''''"'' (A.20c) 

Ply+P2j; = -P;;'^"°'^''^ (A.20d) 

Plz + P2z = sinh(y)sl/2 _ ^branches (^206) 

that is parametrized by the momenta of the branches p^ and p4, by the total momentum 
^branches q£ g^^j ^j^^ branches represented by the blobs in Figure |l^, by the rapidity y and the 
invariant mass s^" of the colliding partons and by the variables -Ei, E2 and p2y The next 
step is to fix the value of the variables Ei, E2 and p2y The variable Ei can be expressed 
as a linear function of E2 : 

El = COsh(y)sl/2 _ ^branches _ ^^^ ^^21) 

Then the mass-shell conditions for the two missing particles define a system of two coupled 
quadratic equations in the variables E2 and p2y In this case, the quartic terms of the two 
equations have the same coefficients, and the system reduces to a linear equation and a 
quadratic equation. There are up to two solutions for E2 and p2y Each solution that is 
physical (i.e., such that E2 > 0, Ei > 0) corresponds to a distinct phase-space point at 
which the Jacobian in Eq. ( [A.19 ) and the integrand must be evaluated. 



A.6 MB F 

The notation for the phase-space variables associated with this constrained sector is given in 



Figure 17, The momenta of the missing particles are denoted by pi and p2, the momenta 
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Figure 17: Notation for the kinematics of MB F. 



of the branches directly connected to these missing particles are denoted by p^ and ^4. 
The Bjorken fractions are denoted by gi and (72- The variables Sjj refer to the invariants 
{pi+Pj)"^- The standard phase-space parametrization associated with this main block reads 

d?pi d?P2 



dqidq2 



(2^)^5^ (Pi, - Pi 



finj 



(A.22) 



' (27r)32£;i {2tt)^2E2 

The four-vectors Pjn and Pa^ refer to the total momenta in the initial and final states, 
respectively. In our procedure we apply a change of variables that leads to the following 
parametrization of the phase-space measure 

1 



Wir'^EiEi 
The jacobian J of this transformation is given by 

E1E2 



dqidq2dsi3ds24 x J- 



(A.23) 



J 



E4,{pizP2yP3x -piyP2zP3x - PlzP2xP3y + PlxP2zP3y + PlyP2xP3z- 

PlxP2yP3z) + E2PlzP3yPAx " EiP2zP3yP4.x " E2PlyP3zPAx + EiP2yP3zP4x 

E2PlzP3xPAy + EiP2zP3xP4y + E2PlxP3zPAy " EiP2xP3zP4y + {E2PlyP3x + 

-Eip2yP3x - E2PlxP3y + EiP2xP3y)PAz + E3{-pizP2yP4x + PlyP2zP4x 

-1 

+PlzP2xP4y - PlxP2zP4y " PlyP2xPAz + PlxP2yP4z 



(A.24) 



If we treat the variables Ei, E2 and p2y as independent parameters, the other compo- 
nents of the momenta pi , p2 of the missing particles can be expressed as the solution of the 
following linear system of five equations. 



{pi+P3f = si3 (A.25a) 

{p2+P4f = S24 (A.25b) 

Pix+P2x = -p'^'^'"''"'' (A.25c) 

Ply+P2y = -P'^''''"''''' (A.25d) 

Piz + P2z = s^/\qi - q2)/2 - p^--'^- (A.25e) 

that is parametrized by the momenta p3 and p4, by the total momentum jf^^'^"^^^ of all the 



branches represented by the blobs in Figure 17, by the Bjorken fractions qi, q2 and by the 



34 



variables Ei, E2 and p2y The next step is to fix the value of the variables Ei, E2 and p2y 
The variable Ei can be expressed as a linear function of E2: 

El = S^/^{qi + q2)/2 - ^branches _ ^^_ ^^26) 

The mass-shell conditions for the two missing particles with momenta pi and p2 define a 
system of two coupled quadratic equations in the variables E2 and p2y In this case, the 
quartic terms of the two equations have the same coefficients, and the system reduces to 
a linear equation and a quadratic equation. There are up to two solutions for E2 and p2y 
Each solution that is physical (i.e., such that £'2 > 0, i^i > 0) corresponds to a distinct 
phase-space point at which the Jacobian in Eq. ( A.24| ) and the integrand must be evaluated. 



B. Phase-space measure associated with the secondary blocks 

B.l SB A 

The notation for the phase-space variables associated with this secondary block is given in 
Figure ll^. The momentum of the missing particle is denoted by pi, the momenta of the 



Figure 18: Notation for the kinematics of SB A. 

three branches connected to the block are denoted by p2, P3 and p^. The variables S12, 
S123 and S1234 refer to the invariants (pi +^2)^, (pi +P2 +^3)^, and {pi +P2 +P3 +P-if-, 
respectively. The standard phase-space parametrization associated with this block reads 

In our procedure we apply a change of variables that leads to the following parametrization 
of the phase-space measure 

/o NQor. dsi2dsi2-idsi2M X J. (B.2) 

[iTryiEi 

The jacobian J of this transformation is given by 

El 

J = -^ EA{pizP2yP3x -PlyP2zP3x - PlzP2xP3y + PlxP2zP3y + PlyP2xP3z 
o 

-PlxP2yP3z) + E2PlzP3yPix " EiP2zP3yPAx - E2PlyP3zPAx + EiP2yP3zPAx 

-E2PlzP3xPiy + EiP2zP3xP4y + E2PlxP3zPAy " EiP2xP3zPiy + {E2PlyP3x 

-EiP2yP3x - E2PlxP3y + EiP2xP3y)PAz + E^{-pizP2yPAx + PlyP2zPAx 

-1 
+PlzP2xPiy - PlxP2zPiy " PlyP2xP4z + PlxP2yP4z) ■ (B.3) 
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If we treat the variable Ei as an independent parameter, the components of the three 
momentum p^ of the missing particle can be expressed as the solution of the following 
linear system of three equations 

{pi+P2f = si2 (B.4a) 

{Pl + P2 + P3f = S123 (B.4b) 

(Pl +P2+P3+ Pif = 51234 (B.4c) 

that is parametrized by the momenta p2 ,P3, Pa of the branches connected to the block, 
and by the variable Ei . The next step is to fix the value of the variable Ei . The mass-shell 
condition for the missing particle with momentum pi defines a quadratic equation in the 
variable Ei. There are up to two solutions for Ei. Each solution that is physical (i.e., 
such that El > 0) corresponds to a distinct phase-space point at which the Jacobian in 
Eq. (|B.3|) and the integrand must be evaluated. 



B.2 SB B 

The notation for the phase-space variables associated with this secondary block is given in 
Figure ^. The momentum of the missing particle is denoted by pi, the momenta of the 
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Figure 19: Notation for the kinematics of SB B. 

two branches connected to the block are denoted by p2 and p^. The variables si2 and S123 
refer to the invariants (^1+^2)^! and (^1+^2+^3)^, respectively. The standard phase-space 
parametrization associated with this block reads 

^'^^ (B.5) 



(27r)32£;i ■ 

In our procedure we apply a change of variables that leads to the following parametrization 
of the phase-space measure 

1 



(27r)32£;i 



!)i(isi2(isi23 X J, (B.6) 



where (pi denotes the azimuthal angle of particle i. The jacobian J of this transformation 
is given by 



T ^1 



COs((/)i - (l)2)E3P2TPlz + COs((/)i - (j)3) E2P3TPIZ + E3P1TP2Z 

-1 
COs{(j)l - (l)3)EiP3TP2z - E2PITP3Z + COs{(pl - (l)2)EiP2TP3z ■ (B.7) 
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If we treat the variable Ei as an independent parameter, the transverse momentum piT 
and the momentum component piz of the missing particle can be expressed as the solution 
of the following linear system of two equations 

{Pi+P2f = si2 (B.8a) 

{Pl + P2 + P3f = Sl23 (B.8b) 

that is parametrized by the momenta p2 and pa of the branches connected to the block, 
by the azimuthal angle 0i and by the variable Ei . The next step is to fix the value of the 
variable Ei . The mass-shell condition for the missing particle with momentum pi defines a 
quadratic equation in the variable Ei . There are up to two solutions for Ei . Each solution 
that is physical (i.e., such that Ei > 0) corresponds to a distinct phase-space point at 
which the Jacobian in Eq. (|B.7D and the integrand must be evaluated. 



B.3 SB C/D 

The notation for the phase-space variables associated with this secondary block is given in 
Figure pQ. The momentum of the missing particle is denoted by pi , the momentum of the 
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Figure 20: Notation for the kinematics of SB C/D. 

branch connected to the block is denoted by p2- The variable si2 refers to the invariant 
{pi +^2)^- The standard phase-space parametrization associated with this block reads 

^'^^ (B.9) 



(27r)32£;i ■ 



In our procedure we apply a change of variables that leads to the following parametrization 
of the phase-space measure 

^ ' hd9idsi2 X J, (B.IO) 



{2tt)^2Ei 

where 9i and <j)i denote the polar and azimuthal angles of the missing particle. The jacobian 
J of this transformation is given by 

1 



J = — sme^ilpil 



\P1\E2 - Eip-^.P2 



(B.ll) 



If we treat the variable Ei as an independent parameter, the momentum modulus \pi\ of 
the missing particle can be expressed as the solution of the following linear equation 

ipi+P2f=Sl2 (B.12) 
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that is parametrized by the momentum p2 of the branch connected to the block, by the 
polar and azimuthal angles 9i, (pi, and by the variable Ei. The next step is to fix the 
value of the variable Ei . The mass-shell condition for the missing particle with momentum 
pi defines a quadratic equation in the variable Ei. There are up to two solutions for Ei. 
Each solution that is physical {i.e., such that Ei > 0) corresponds to a distinct phase-space 
point at which the Jacobian in Eq. (B.ll) and the integrand must be evaluated. 



B.4 SB E 

The notation for the phase-space variables associated with this secondary block is given in 
Figure ^. The momenta of the visible particles are denoted by pi and p2, the momentum 
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Figure 21: Notation of the kinematics for SB E. 

of the branch connected to the block is denoted by p^. The variables S12 and S123 refer to 
the invariants (pi +^2)^, and (pi + P2 +^3)^- The standard phase-space parametrization 
associated with this block reads 



d^pi (fp2 



(B.13) 



(27r)32Ei (27r)32^2 ' 

In our procedure we apply a change of variables that leads to the following parametrization 
of the phase-space measure 

1 



-d9id(j)id92d(f>2dsi2dsi23 x J, 



(B.14) 



{27r)^4:EiE2 

where 6i and (pi denote the polar and azimuthal angles of particle i. The jacobian J of this 
transformation is given by 



j^2 

J = —r\Pi\ sin0isin02 



\Pi\\P2\/Ei - \P2\fi2){E3 - \P3\f23) 



-{E^Ml/Ei - \p^\fi3){Ei - fi2\Pi\ 



(B.15) 



where fij stands for Pi-Pj/\Pi\\Pj\- The values for the momenta \pi\ and IP2I can be 
obtained by solving the following linear system of equations 



{Pl+P2f 
{Pl +P2+P3f 



S12, 

Sl23- 



(B.16a) 
(B.16b) 



By subtracting Eq. ( |B.16b ) from Eq. ( |B.16a ), we obtain an expression for Ei that is a first 
order polynomial in \pi\ and \P2\- Inserting this expression into Eq. ( [B.16a| ) and into the 
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equation defining the mass shell condition for the particle of momentum pi, we obtain a 
system of two quadratic equations in |p^| and IP2I parametrized by the momentum ^3, by 
the invariants S12, S123 and by the angles 61, 62, (pi and 4>2- This system can be solved 
analytically. There are up to four solutions for the modulus |p^| and \P2\- Each solution 
that is physical (i.e., such that \pi\ > and IP2I > 0) corresponds to a distinct phase-space 
point at which the Jacobian in Eq. ( B.15| ) and the integrand must be evaluated. 
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